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ABSTRACT 


The dynamic analysis of large, complex structural systems is 
computationally intensive and therefore prohibits the use of optimization 
procedures, which are both iterative and complex with respect to variable 
search patterns. The solution to this problem is through the use of time and 
frequency synthesis techniques. They provide a means of rapidly recalculating 
a system’s changed response due to structural modifications, as dictated by the 
optimization procedure. The efficiency is gained through the fact that the 
synthesis methods are independent of model size, in that only those model 
degrees of freedom where changes are made are required in the analysis. 
Furthermore, these methods are exact in their formulation, including the 
treatment of non-proportional damping. These structural synthesis 
techniques are developed in the context of optimal design of shock and 
vibration isolation systems. Their utility and value is demonstrated in the 
optimal design of an isolation system for a 109 dof non-proportionally 
damped structural system. In the course of the optimization, the synthesis 
techniques make possible 80 transient, frequency response, and static analyses 
in 2 hours and 39 minutes (desktop computer), while yielding an isolation 
design which satisfies all design constraints. 
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I. INTRODUCTION 


With the development of the finite element (FE) method, and the 
advancement of computer technology, it is now possible to design and 
analyze complex engineering systems. Through FE techniques, a detailed 
mathematical model of a complex structural dynamic system may be 
developed, and the static and dynamic responses determined. The 
‘traditional’ procedure for conducting a finite element analysis (FEA) of a 
structure is to first assemble the system matrices. These system matrices are 
the mass, stiffness, and damping matrices, which constitute the FE model of 
the complete structure. The next step would be to determine the system 
responses using various solution techniques. This is the most 
computationally demanding phase of the FEA process. The results are then 
processed and interpreted. As a result of this analysis, the engineer may wish 
to change some aspect of the design and perform the analysis again. Even 
more useful would be the use of some type of optimization scheme to find 
the optimum design change while concurrently performing the FEA. 
However, if the ‘traditional’ method of FEA is used, every time a design 
variable is changed, the affected system matrices must be reassembled and the 
entire solution phase must be reaccomplished. Depending on the complexity 
of the system, and the number of different design parameters which may be 


changed, this route would be computationally impractical. As a result, the 


designer may only be able to iterate through the design process a few times, 
and be left with a design which is less than optimum. 

Therefore, more efficient techniques for calculating a system's 
responses after modifications have been made, must be introduced. One such 
technique is the use of synthesis procedures to obtain the modified system 
responses. These synthesis procedures involve the use of the structure’s 
baseline (pre-modification) responses, the modifications made to the 
structure’s mass, stiffness, and damping matrices, and the equations which 
link the two to obtain the modified/synthesized responses. The advantage of 
this is that the system matrix assembly and solution phase of the FEA process 
must only be accomplished once. A smaller set of change matrices are 
assembled, and along with the presynthesized responses and computationally 
efficient synthesis equations, the synthesized responses are obtained. 

The two types of sythesis techniques to be featured in this thesis are 
frequency domain synthesis, and time domain synthesis. Frequency domain 
synthesis makes use of the baseline frequency response functions (FRFs) to 
calculate the new FRFs after the structure has been modified. Time domain 
synthesis uses the impulse response functions (IRFs) of the baseline structure, 
the coupling forces from the modification, and the convolution integral. The 
combination of these produces a nonstandard nonhomogeneous Volterra 
integrodifferential equation (VIDE) of the second kind which is solved in 


order to calculate transient responses of the modified structure. The use of 


these techniques greatly reduce the computer computation time and allow for 


the application of optimization techniques to complex engineering systems. 





II. FINITE ELEMENT FORMULATIONS 


As mentioned previously, the initial step in the FEA process is to create 
a finite element model of the structure to be analyzed. In this study, different 
finite element types are used to ensure the validity and universal applicability 
of the synthesis techniques, and to assist in the synthesis formulations. For 
simplicity, linear lumped mass-spring systems were used. Its formulation is 
based on the use of Newton’s laws to derive the equations of motion for each 
particle of mass [Ref. 1]. These systems allowed for the initial general 
formulation of the synthesis techniques and the building of a checking system 
to ensure its accuracy. In order to illustrate more realistic systems, systems 
formed from beam elements and plate elements were used. The beam 
element formulation is based on Euler-Bernoulli beam bending, and the use 
of Hermitian shape functions [Refs. 2 & 3]. The plate element formulation is 
based on shear deformation which includes both the transverse shear energy 
and the bending energy [Ref. 3]. One very important anomally about the plate 
shear deformation formulation is that, the unconstrained stiffness matrix (K) 
is rank deficient by the number of degrees of freedom per node, plus one. 
When solving the eigenvalue problem to obtain the system’s natural 
frequencies (@,) and mode shapes, one additional “spurious” rigid body mode 
is obtained which disrupts all modal calculations. Therefore, to avoid this 


problem, all analyses done using the plate were done on a constrained to 


ground system to eliminate all rigid body modes. The constraints were then 
subsequently synthesized out during the analysis in order to obtain the truly 


desired responses. 


IJ. FREQUENCY DOMAIN STRUCTURAL SYNTHESIS 


The following background information on frequency domain 
structural synthesis is taken from references [4 & 5]. This portion of the thesis 
will outline the methodologies and formulations of the frequency synthesis 
technique and the classical techniques used to check its accuracy. 

Frequency domain techniques were demonstrated as early as 1939 by G. 
Kron in his book on the tensor analysis of electrical networks. Since then, 
numerous formulations and application techniques have been used and 
documented. The advantages in computational time, and the flexibility of 
this technique, have also been well documented [Ref. 6]. The solutions 
obtained through this method are exact with no approximations. 

As mentioned previously, frequency domain structural synthesis is a 
methodology whereby changes can be made to a given structure, and its new 
frequency response function (FRF) can be formulated through the use of a 
single synthesis equation. This is quite different from the classical method of 
evaluating a structural change through the reconstruction of the basic 
elements which define the structure, and then using an impedance inversion 
or modal superposition technique on the new structure. The frequency 
synthesis technique may be divided into two major classes, coupling and 
modification. Coupling is the joining of a totally independent uncoupled 


structure to the given structure. Modification is the addition of redundant 


load paths in the given structure. Each of the two classes of synthesis may be 
either direct or indirect. For the indirect case, there is the existence of an 
interconnection impedance element between substructures in a coupling 
operation, and degrees of freedom in a modification operation. In the direct 
case, the coupling and modification is done without the impedance element 
existing, therefore making the modification synonymous to applying a 
constraint. 

This thesis will focus on the indirect modifications to a given structure. 
It will allow for the inclusion of a spring stiffener or a viscous damper 
between two points, a lumped mass on the structure, and implementing a 
base excitation to the structure. A generalized computer program will be 
presented, which will allow for any of the above mentioned changes to be 
accomplished on some baseline structure. The program then returns the FRFs 
for all dofs where changes have occurred and any other dofs that may be of 


interest to the user. 


A. GENERALIZED FREQUENCY RESPONSE 

The initial step in obtaining the new synthesized FRFs for the modified 
structure, is to have the FRFs for the baseline structure. The following 
equations are for a system which is excited by a force on the structure or by a 
base excitation. These formulations are also used to verify the accuracy of the 


synthesis after a modification has been made. 


Consider the following two degree of freedom system: 





x(t) x(t) 


Figure 3.1 Two Degree of Freedom Mass-Spring-Damper System 


where: 
m, = mass of a plate 
m,. = mass of computer equipment 
k, = isolator stiffness between plate and ground 
k. = isolator stiffness between plate and computer equipment 
x,(t) = plate displacement as a function of time 
x(t) = computer displacement as a function of time 


F(t) = excitation force 


Applying Newton’s 2nd law, the equations of motion for this system may be 


written as: 


0 ia A * +, She rn r +k. rit . {0| G1) 
O m, j|x, ce Case ake k, |x, 0 
Generalizing this formulation to an ndof (number of degrees of freedom) 
system where the external force may occur at any dof, and using a more 
compact matrix and vector notation, the 2nd order system of linear equations 
for an ndof structure can be written as: 
[M]ix} +[C]ix} + [K]ix} = {F} (3.2) 
Assuming a harmonic forcing function and consequently a harmonic steady- 
state response: 
{F()} ={F }e™ {x(t)} ={X}e™ (323) 
where: 
{F\ & {X} are the amplitudes of the harmonic forcing function and 
response. 
Taking the first and second derivatives of the response {x(t)}, and substituting 
into equation (3.1) and simplifying: 
[K+ jQC-Q?M]{X} ={F} (3.4) 
[Z(Q){X} ={F} (3.5) 
[Z(Q2)] is known as the impedance matrix. Multiplying both sides of the 


equation by [Z(Q)]’ and denoting this as frequency response function (FRF) 


10 


matrix [H(Q)], equation (3.5) becomes: 


{X}=[H(Q)]{F} (3.6) 
where: 
Ay Fi, cee nor 
| ito (3.7) 
Baier tmetndorx2e nace endos 
and in elemental form 
H, -? (3.8) 


The FRF matrix [H(Q)] is both complex valued, and frequency dependent. The 
advantage of this formulation in this work is that it can be used to check the 
frequency synthesis method when damper modifications are made to the 
structure. These damper modifications represent non-proportional damping, 
which is not easily handled in modal coordinates. 

The modal coordinate formulation of the FRF matrix begins from the 
same general equations of motion presented in equation (3.2), and the 
assumptions of harmonic excitation and response in equation (3.3). It also 
includes the assumption of a harmonic modal response: 

19) } = {7 her (3.9) 


and the linear modal transformation: 


{x(t)} =[®]{q@} (3.10) 


I] 


where [®] is the assembles matrix of mass normalized mode shapes. By 


applying equations (3.3), (3.9), & (3.10), to equation (3.2) and simplifying: 
|-Q?[M][®] + jOlC]]+[K]®]|{7} = {F } (3.11) 


Premultiplying equation (3.11) by [®]’: 


—< 1 44+9jQ) 82a +| @, {q}={} (3.12) 


where it is recognized that: 


[®]’ [M] [] = [1] (3.13) 
[D}*[C] [©] = [2a (3.14) 
[®]’ [K] [©] = [ol (3.15) 
[D]’ {F} = {3} (3.16) 


¢ is the modal damping factor, @ is the natural frequency, and the subscript r 
represents each mode. Rewriting equation (3.12) and solving for the modal 


displacement yields: 


iq}= {g} (3.17) 


w? —Q? + 526.0, 


Using equations (3.10) & (3.16), equation (3.17) is transformed back to physical 


coordinates: 


I 


{X}=[®] | 


or O+peon fi oe 
From equation (3.6): 
fH(Q)} = [0] : [o}" (3.19) 


w.” -Q? + 526.0. 


and any element of the FRF matrix may be represented as: 


a 9,9) 


Although this modal coordinate formulation does not handle non- 
proportional damping, it does provide for the ability to sum over a subset of 
the complete modes, rather than summing over all modes for the FRF. The 
number of modes necessary to correctly capture the response is dependent on 
the frequency range of interest. The lower the frequency range of interest, the 
less number of modes necessary. The obvious advantage of this procedure is 
that for large dof systems, this summation may be truncated, and 
computational time saved. However, the question of how many modes are 
enough must be considered and truncation criteria must be established. The 
modal formulation will be the one used for the synthesis techniques. 
Recognizing the modal truncation issue, this thesis uses all modes in its FRF 


calculations. 
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Now consider the following two degree of freedom system: 





y(t) x, (t) x(t) 


Figure 3.2 Base Excited Two Degree of Freedom Mass-Spring-Damper System 


where: 
k, = isolator stiffness between plate and base 
y(t) = base excitation displacement 


all other variables are the same as in the system illustrated in Figure 
Sal) 


Applying Newton’s 2nd law, the equations of motion for this system may be 
written as: 
m, 90 |x, Cp Co SC. IX, |__| Kpteeeercen | |S |  Commaelli k, 0 
a mile(*] - oc lxft] ke ok xf lo of!t1o of 
(3.21) 
Generalizing this formulation to an ndof system where the base excitation 


may occur at any dof, and using a more compact matrix and vector notation, 


the 2nd order system of linear equations for an ndof structure can be written 
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as: 


[M]iX} +[C]HX} +[K]x} =[C. Hiv} +[K. hiy} (3.22) 
Following the same procedure as with the FRF formulation for a force 


excitation on the structure, and assuming a harmonic base excitation 


({y(t)} ={Y $e), the following result is obtained: 


{x}=[K+jQC-Q’M] [K, +72C,]{¥} (3.23) 
Since the base excitation iY} is the same at all dofs, Y can be factored out and 


the result is: 


{x} =[K+jQC-Q°M] [K, +72C, |{1}Y (3.24) 
and 
{H}=[K+jQC-Q°M] [K, +j2C,}{1} (3.25) 


In this instance the collection of FRFs is not a matrix but rather a vector 


where in elemental notation: 


H, = 


1 


(3.26) 


|| ><! 


In some instances, the base excitation and desired response may not both be 
displacements. One different combination could be an interest in obtaining 
the output displacement response due to a given base acceleration. From the 
assumption of a harmonic base excitation and response: 


Ly(t)} = {¥ be™ {x(t)} = {Xbe* (3.27a&b) 
Ly} = {¥ te = -071Y te™ {x(t)} = {x}e™ =-O7{X}e™ — (3.28a&b) 


ls 


» {Y= (a |{* {xX} = steaite (3.29a&b) 


The following table uses the above relationships between acceleration and 


displacement to illustrate possible combinations and how to handle each. 





Table 3.1 Response and Base Excitation Relationships 
The modal formulation for base excitation is not needed since the 
synthesis method for a base excitation or force excitation, uses the previous 


modal FRF formulation exclusively. 


B. FREQUENCY DOMAIN SYNTHESIS FORMULATION 


Now consider the following general ndof (number of degrees of 
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freedom) structure, to which structural modifications are to be made: 


~~ 


Figure 3.3 General NDOF Structural System 


The letters 1, c, and b represent the physical coordinate system for the 
structure where: 

1=iset - the set of internal dofs where no changes occur, but there is an 

interest in knowing FRF information about these dof. 
c=cset - the set of dofs where changes have occurred. 
b= bset - the set of dofs where the structure is indirectly connected to a 
base excitation. 

Using equation (3.6) and the previously described coordinate system for 
rearrangement and partitioning, the structural system is described in the 


frequency domain at each frequency as: 


x? =| [Ha] ‘TA.] ; [Ha] hf (3.30) 


Xp [Hi [H..] : [Hye | |lé 
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Before modification, the force vector refers to externally applied forces 
to the structure. However, following modification, there exists coupling 
forces between the cset and bset dofs. By definition, the iset dofs do not 
experience these coupling forces. Therefore, the partitioned forces may be 


rewritten as: 


f =o (3.31a) 
f =f +f? =f —AZ x’ (3.31b) 
f, = fo +f? =f — AZ, (x? -y) (3.31) 


where: 

f* - the coupling forces due to structure modification. 

[AZ] - the impedance matrix which is equal to [AK + jQAC + Q’?AM]. 

x" - generalized (synthesized) responses after modification. 
Substituting equations (3.31) into equation (3.30), and recognizing that the 
presynthesized responses {x, x. x,}' are those due to the external forces only, 
the following result is obtained: 

x:] (x, [Hi] [Hi] [AZ] 0 {x [H;.] [Hs] [AZ.] 0 fo 
xf =4xe¢-|[H] [Ha] 0 [az it \. [H..] [Ha] ce zat 
J 06) [Hoc] [ue] mee THe] [Hee] 
(3.32a) 
Denoting two additional coordinate sets ‘e’ and ‘z’ where: 
e =eset =1UcUb 


z=zset=cUb 
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and [AZ] is the total impedance matrix. Equation (3.32a) can be rewritten as: 

1Xe } = {Xe} y [FH., [AZ] x, } + [H., |[4Z, Ky} (3.32b) 
From equation (3.32b) it can be seen that there are two unknown synthesized 
responses in this one equation. Another independent equation can be 
generated from equation (3.32a) without introducing any additional 
unknowns. By observing the bottom two rows of equation (3.32a), the 
following relationship is derived: 

{x.} = {x.}-[H[AzZ]x.} +[H. [42 Hy} (3.33) 
Solving for {x,}" and substituting that expression into equation (3.32b) yields 


the general expression for the responses of the synthesized structure. 


{x.} ={x.}-[H.[Az|([I+H,,42]"{x,}+[1+H_,42] [H,]4Z, ]{y}) +[Ha [42 Hy} 
(3.34) 
From the general equation for frequency response, it is recognized that: 
{x.}=[H.ef"} {x,}=[H,}&"} (3.35a&b) 
Equation (34) is therefore rewritten as: 
{x.} =[H.]{f"}- [He [az ]((l+ 42) [H,.]{f"} + [+ Haz) [H,. 42, ]{y}) + [Ha [420 ]{y} 
(3.36) 
As can be seen from equation (3.36), the synthesized response is now in terms 
of the known FRFs of the presynthesized structure, known impedance due to 


modifications, and known force and base excitations. 
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Now that the complete general frequency response synthesis equation 
has been derived, consider the case where there is no base excitation applied 
to the structure. By letting {y} = 0, the general synthesis equation becomes: 

{x,}' = ((H..]-[H..[4Z][1+ H,,2] [H,,]}{"} (3.37) 
Since no base excitation exists, the impedance between the structure and base 
(AZ), is reduced to a ’c’ dof to ground change. Therefore the bset dofs are 
really additions to the cset: z=cUb =cUc =cand AZ, = AZ.. Also 
recognizing that {x,} = Eee the synthesized FRF matrix at each 
frequency can be written as: 
fH.) =[H.J- [Haz J+ H,.02,}"[H.] (3.38) 
Now consider the case where there are no external forces applied to the 


structure. By letting ee = 0, the general synthesis equation becomes: 


(x.}" = ((H.,[4Z,]-[H.. [az ]0+H, a2] 1H, [4z,]){y} (3.39) 
Factoring out the base excitation magnitude Y and applying equation (3.6) yet 
again: | 

(H,}' = ((HJ[4Z,]-[H..[42]0+ 8, 42] [H,J4Z,)}} 40) 
In elemental notation, the FRF vector represents the following at each 
frequency: 


fH} ={x, % (3.41) 
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Equations (3.38) and (3.40) show how the new synthesized FRFs are 
obtained from the presynthesized FRFs and the impedance caused by 
structural modifications. These synthesis equations are exact, and there are no 
limitations to the modification values. The change may even be negative to 
represent the removal of mass or removal of an interconnection element 


from the structure. 


eo FREQUENCY DOMAIN SYNTHESIS COMPUTER CODE 

The computer language used for the frequency domain synthesis and 
all other computer coding in this thesis is MATLAB V.4.2c.1. The goal in the 
formulation of the frequency domain synthesis computer programs is to 
perform the synthesis for an indirect modification operation to a given 
structure. The program allows for the inclusion of a spring stiffener or a 
viscous damper between two dofs, a lumped mass addition on the structure, 
and implementing a base excitation to the structure. The program then 
returns the FRFs for all dofs where changes have occurred and any other dofs 
that may be of interest to the user. Appendix A shows the computer codes 
used to perform the frequency domain synthesis, and perform a comparative 
analysis of the synthesis versus classical methods. The term classical is 
synonymous to the traditional method of reformulating the elemental 


matrices following a modification. 
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The programs presented in Appendix A can be segregated into three 
categories. The first is the calculation of the synthesized FRFs due to an 
external force excitation on the structure, and their comparison to the 
classically calculated FRFs. The second is the calculation of the synthesized 
FRFs due to a base excitation on the structure, and their comparison to the 
classically calculated FRFs. Both of these programs use the impedance 
inversion method to calculate the presynthesized FRFs, and the classical FRFs 
following modification. The third category is the calculation of the 
synthesized FRFs where the different aspects of the synthesis technique have 
been modularized into MATLAB functions. This modularization assists the 
process in running more efficiently, and allows for the universal use of the 
frequency synthesis technique. This is crucial since the ultimate goal of the 
technique is to be used efficiently in another process (i.e. optimization). Since 
efficiency is of concern, the modal method of calculating the presynthesized 
FRFs is used. 


The following is a diagram of the flowpath of the modularized 


z2 


frequency synthesis programs: 


Load FRFs 

change.m 

Load Structure Calculate FRFs Modifications Form Change Perform Synthesis FRF Selection 
Matrices 

change.m modal.m change.m frfsynth.m frf_sift.m 


frfmodal.m deltam 





Figure 3.4 Flowpath of the Modularized Frequency Synthesis Programs 


A more detailed description of each function is contained within Appendix 


A. 


D. ILLUSTRATIVE EXAMPLES 

The purpose of this section is to illustrate the use of the frequency 
domain synthesis technique and validate its accuracy. The computational 
efficiency of this technique has already been well documented [Ref. 6], 
therefore a time comparison is not done. The three structural systems 
previously described, mass-spring-damper, beam and plate, are used to 


accomplish this goal. 


ZS 


i Mass-Spring-Damper System 


Consider the following mass-spring-damper system: 


k,, 
pe k = 1000 lb/in 
k ! ; k,, = 250 Ib/in 
FN FR II m = 0.259 Ib-s’/in 
Z fa ali F(t) Cy) = 2 Ib-s/in 
l C36 
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Figure 3.5 Mass-Spring-Damper System Experiencing Force Excitation 


The solid elements represent the original structure and the dashed elements 
represent modifications. The subscript ‘i’ in the excitation force represents the 
dof where the excitation is applied. Although not shown in the figure, the 
original structure is proportionally damped using [C] = a[K]. The 
modifications are arbitrarily selected and the addition of the damper element 
represents nonproportional damping. The following is the resultant 
synthesized and classically determined FRF H,.(Q) which represents the 


response magnitude at dof 3 due to a force at dof 2. 
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Figure 3.6 Force Excitation Mass-Spring-Damper System FRF 


From the graphs produced, it can been seen that the plots are identical, 
proving that the synthesis technique is exact, and that the computer coding 


for the arbitrary changes is correct. These results were produced using the 


fsynstr.m program. 
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Zz: Beam System 


Now consider the following beam system: 


k, = 1000 lb/in 
C139 = 20 Ib-s/in 
m,, = .5 lb-s’/in 





Figure 3.7 Beam System Experiencing Base Excitation 


where the beam parameters are : 

Length: L=5 ft Width: w = 3 in Height: h = 4in 

Young’s Modulus: E = 10e6 psi 

density: p = 2.53e-4 lbf-sec42/in*4 
The free-free beam which is discretized into 10 beam elements, 11 nodes, and 
22 dofs, represents the original structure. The modifications are represented 
by the dashed lines. Just as in the previous example, proportional damping is 
used to form the original [C] matrix and all changes to the beam are arbitrary. 
The following is the resultant synthesized and classically determined FRF 
H,,(Q) which represents the response magnitude at dof 13 due to the base 


excitations. 
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Figure 3.8 Base Excitation Beam System FRF 


Just as before, it can be seen that the two plots are identical. These results were 


produced using the fsynbase.m program. 


2 Mass-Plate System 

Up until now the structures use to demonstrate the frequency synthesis 
technique were relatively simple with relatively small number of dofs. The 
true test of the technique and the computer code is in their applicability to 


relatively large structures. Therefore, consider the following computer-plate 
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system: 






Ak, = 1e4 Ib/in 
Ac, = 1 lb-s/in 
Ak, = 500 Ib/in 
Ac, = 1 Ib-s/in 


c, | 


mf y(t) Ji a 


Figure 3.9 Plate-Mass System Experiencing Base Excitation 


where the computer and plate parameters are : 
Plate: Width: w =5 ft Depth: d =5 ft Thickness: t = 1 in 
Young’s Modulus: E = 30e6 psi Poisson’s Ratio: v = 0.3 
Density: p = 7.35e-4 Ibf-sec42/in”4 
Computer: Weight: W = 100 lbs 
In this case, the original structure is represented by the free-free plate 
with the attached computer located off center and above the plate. The plate is 
discretized into 25 plate elements, 36 nodes, and 108 dofs. The computer is 
modeled as a single lumped mass in the translational direction only, and 
increases the total number of nodes and dofs to 37 and 109 respectively. The 


modifications are represented by the dashed lines and are comprised of 
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changes to the isolators between the base and the plate. Just as in the previous 
example, proportional damping is used to form the original [C] matrix, and 
any changes to the system are arbitrary. However, since the manipulation of 
this system in an optimization routine is the ultimate goal of this research, 
the isolator elements which will act as optimization variables are chosen for 
modification. Therefore, modifications were also performed on the isolators 
between the plate and the computer. The modularized program which uses 
the modal method of calculating the original FRF will be used in this 
instance. However, the classical method will also be used as a check system. 
The following is the resultant synthesized and classically determined FRF 
H99(2) which represents the response magnitude at dof 109 (the computer) 


due to the base excitations. 
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Figure 3.10 Free-Free Plate-Mass System FRF (Modal Synthesis) 


From the plots it is easily observed that the synthesis method and the 
classical method do not yield the same results. After successfully performing 
countless other frequency synthesis analysis on all three structures, and 
successfully re-running this same analysis where the original FRF is obtained 


using impedance inversion methods vice modal methods (Fig. 3.11 shows), 
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Figure 3.11. Free-Free Plate-Mass System FRF (Impedance Synthesis) 


it is concluded that the problem lies in the use of the modal method of 
calculating the original FRF. After observing the mode shapes and natural 
frequencies of the original free-free plate-mass system, the existence of an 
additional spurious rigid body mode is ered. The existence of this rigid 
body mode vice a flexible mode makes the FRF calculations incorrect. After 
unsuccessful attempts at trying to replace the rigid body modes naturally 
generated by this finite element formulation, with rigid body modes created 
using the Graham-Schmidt method [Ref. 2], the decision to constrain the 
plate-mass and then remove the constraints using synthesis is reached. What 
this accomplishes is to eliminate the rigid body modes entirely and allows for 


the proper calculation of the original FRF. However, there is a computational 


Sit 


price to pay with the addition of four more load paths to contend with in the 
calculations. This is however, a very minute price to pay when compared to 
using the impedance inversion method. 

The following is the resultant modally-synthesized and classically 
determined FRF H,,.(Q), when the original plate-mass structure is exactly as 
before, but with 1000 lb/in spring-to-ground elements located at all four 
corners. In the process of the analysis, these elements are synthesized out of 


the structure. Other modification values are still the same as before. 
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Figure 3.12 Corners-to-Ground Plate-Mass System FRF (Modal Synthesis) 


It can be seen from these plots that eliminating the rigid body modes in the 


original structure, and then synthesizing the constraint elements out is an 
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accurate means of avoiding the problem with the rigid body modes. Another 
verification of the ability to synthesize out elements is that fact that these 
plots also match the plots in Figure 3.11 (free-free plate-mass system); as well 
they should. This will therefore be the method used in order to handle the 


dynamic analysis of the plate-mass system. 


oo 





IV. STATIC DISPLACEMENT SYNTHESIS 


When using optimization procedures for the design of a system, there 
are usually various aspects of the design problem which are at odds with one 
another. For example, in the design of a support beam, the aspect to be 
optimized may be the weight of the beam. One method to reduce weight 
would be to reduce the beam’s cross sectional area. However, since this is a 
support beam and will therefore be carrying some sort of load, a reduction in 
cross section results in an increase in stress. Since there are material 
limitations on the amount of stress the beam can experience, the reduction in 
cross sectional area is constrained by the stress limit. 

For a shock and vibration isolation system, one aspect of the design 
which should be considered is the static sag or displacement of the system. 
The static sag could be used in the optimization problem as a constraint 
which limits the amount the isolator stiffness may be modified. Since the 
solution for the static response of a structure involves the inverse of the 
stiffness matrix ({K]’), for large structures, it would not be computational 
efficient to perform this operation every time the optimization variables are 
changed. Therefore, some sort of synthesis technique must be employed to 
alleviate the need to calculate [K]’ directly. This portion of the thesis explores 
the methodologies and formulations for the use of frequency domain 


synthesis in the calculation of static displacement. From henceforth, this 
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technique will be referred to as static displacement synthesis. Classical 
techniques such as Guyan model reduction are used to check the accuracy of 


these formulations. 


A. GENERALIZED STATIC DISPLACEMENT 

From Newton’s 2nd law, the static equations for an ndof system is: 

(F} = [K] 6 (4.1) 
By defining the force vector as the weight of the structure, and assuming that 
there are no external forces on the structure, the displacement vector {x} 
represents the static displacements of the structure due to its own weight. 
Therefore, the static displacements of the structure due to its own weight 
equals: 

(Jae = IKT'{F} = [KI"IMI{g) (4.2) 
where {g} represents a vector where the gravitational constant appears at all 
dofs that are affected by gravity (i.e. vertical translational dofs). 

Equation (4.2) returns the static displacements at all dofs of the 
structure. The size of the stiffness matrix of which the inverse is taken is ndof 
by ndof. Since the repetitive calculation of [K]’ is unsatisfactory due to the 
possibility of its large size, and a relatively small subset of the static 
displacements of the structure are desired, the use of model reduction and 


Static condensation schemes are investigated. 
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B. GUYAN REDUCTION / STATIC CONDENSATION 

The following concepts and equations for model reduction are taken 
from references [2 & 7]. The Guyan reduction defines a transformation matrix 
[T] which can operate on a subset of the structure’s displacements and return 


the full ndof set of displacements for the structure. 


1 =I) (4.3) 


where: 


[T] a (4.4) 


a =aset - those dofs arbitrarily selected as the set of active dofs. 
o=oset - the remainder of the dofs omitted from the aset. 
Substituting equation (4.3) into equation (4.1), and premultiplying both sides 


of the equation by [T]’ yields: 


{F} = [K]{x, } (4.5) 

where: 
[K] =[T]' [KT] (4.6a) 
(F} =[T]' {F} (4.7a) 


Therefore, the static displacements of the desired dofs equal: 


{x,}=[KI{F} (4.8) 
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It can be seen from equation (4.8) that the main cost of obtaining the desired 


static responses of the system is now the inverse of the reduced stiffness 


matrix {K] which has size of aset by aset vice the inverse of the full stiffness 
matrix which has size ndof by ndof. There is a definite improvement in 
computational efficiency, but the issue of forming the transformation matrix 
[T], which contains the inverse of the sub-matrix [K,,] still exists. For a large 
structure, the size of the aset will usually be substantially less than the size of 
the oset, thereby making [K,,]’ a time demanding operation. In an 
optimization procedure, [K,,]' can be performed prior to the iterations begin, 
and passed into the iteration portion of the optimization program. Even 
though [K,,]"’ must only be performed once, the formation of [T] is still 
undesirable due to the possibility that FRF data is readily accessible and the [K] 
matrix is unavailable. Therefore, a method of obtaining the reduced stiffness 


matrix and reduced force vector from FRF data is desired. 


« STATIC DISPLACEMENT SYNTHESIS FORMULATION 

During the derivation of the frequency domain synthesis, the eset 
coordinate set was defined as iUcUb. This means that the eset includes the set 
of all dofs that the user is interested in, dofs where modifications occur, and 
dofs where a base excitation is attached. In other words, when looking at the 


structure as a whole, eset is synonymous to the set of active dofs or aset. 
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1. Formulation of Reduced Stiffness Matrix [K*] 
Using the aset/oset coordinate system previously described but with 
the eset equaling aset, the expanded version of equation (4.6a) is written as: 
[K] =|K.. -K..K,, ‘K..| (4.6b) 
Now using the same eset/oset coordinate system and the fact that [H] = [Z]", 


the following relationship is written: 


Fe ee ete] -( ea (49) 


Fel! HZ) 1Z1) [0 1 
Carrying out the matrix multiplication of equation (4.9), the first row yields 
the following two relationships: 
H,.Z.. +H,,Z,. =I H..Z,, + H..Z.. =0 (4.10aé&b) 

Solving for H,, in equation (4.10b) and substituting it into equation (4.10a) 
yields: 

HAH 2Z2 026 =| (4.11) 
Solving for H,, yields: 

[H..] =[Zee— ZeoZoo Zoe] (4.12) 
Recognizing that Z = K +jQC - QM, and applying the static condition that 


Q=0: 


[H..(0)] =[K..—KuKoo Kee]. (4.13) 
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From equations (4.6b) & (4.13), it is determined that: 

[K]=[H,.]" (4.14a) 
Therefore, since the reduced stiffness matrix can be determined from an FRF 
matrix, when a structure is modified, the new reduced stiffness matrix can be 


determined from a synthesized FRF matrix. 


[K']=[H.(0)| (4.14b) 


Ds Formulation of Reduced Force Vector {F’} 
Recognizing that {F} = [M]{1}g the expanded version of equation (4.7a) 
in eset/oset partitioning is: 
{F} = IM. —K..K,. "M. | M.,. — K.oKoo Moo {1g (4.7b) 


Another very important concept of Guyan reduction is the reduced mass 


matrix. 
[M] =[T]' [MIT] (4.15a) 
=(M,.-K,.K,. Me. : M,, - Keke Mel ac Si | (4.15b) 
By comparing equations (4.7b) and (4.15b), it can be seen that: 
{F} = [M]{1.}¢ (4.16) 


if and only if 


Laci, ad=ntd= {| an 
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where {1,} is an eset length vector of all ones, and {1,} is an oset length vector 
of ones and zeros, with ones at the translation dofs and zeros at the rotations. 
In effect, what equation (4.17) is stating is that for every row of the 
transformation matrix, the sum of its columns is exactly one or zero. In 
general, this occurrence is not true. It is very dependent on which dofs are 
chosen as the eset. Therefore, it must be determined which dofs are allowable 
choices for the eset in order for this formulation to work. 

Consider the general stiffness matrix of a restrained structure: 

K* = KY + K® (4.18a) 
where K” represents the stiffness matrix of the unconstrained structure and 
K® represents the grounded stiffness matrix. Equation (4.18a ) can be rewritten 
in eset/oset partitioning format as: 


IK®]= eet fe [Ke ot) 


[Ke] | [K. 


vos em ee hee me 


0 
7a (4.18b) 
| 1 


The zeros in the grounded matrix are due to the fact that grounds produce no 
off diagonal information. Now a displacement vector is chosen which 
satisfies the static equations of equilibrium where only forces appear in the 
active or eset dofs, and which produces no strain energy internal to the 


restrained structure. 


ree ete tte pete tot 
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From equation (4.19) it can be seen that: 


He ta 


The vector which satisfies equation (4.19) is already chosen by default. The 


t ¥e 
vector y 


oO 


(eee er 


1 
equals the previously defined 1 


oO 


Now consider the following portion of equation (4.19): 


(4.20a) 


‘ Therefore, {¥,} = {1,} and 


Ie, ae F (4.19) 
| ra ef [0 
Repartitioning this equation into translational and rotational dofs: 
[Ke], 9 = 0 OE) LEP}, 
G 
0 [Koh  ? 0 [t¥ele| _ HE}, (4.21) 
0 OP Ke Oe ea 0 
| 0 0 0 [KEL Per 0 
where: 
(eee {1.}, 
Vie ie ae ee ie (4.22) 
(Yoir}  |{o}s 
(Foie) | {0} 
Case I - Taking the third row of equation (4.21) yields: 
[Koo] {obs = [Ko], toby = {0} (4.23) 


Since [KS] is a diagonal matrix, the only way equation (4.22) is true is if 


Ke |, =[0]. Therefore, this proves that there cannot be any translational 


42 


grounds in the oset. Consequently, this means that all translational grounds 
must be in the eset. 


Case II - Taking the fourth row of equation (4.21) yields: 
KS], {¥.}, = [KS], {0} = {0} 4.24) 
Equation (4.24) is always true regardless of [Sak Therefore, this shows that 


there are no restrictions on rotational grounds being in the oset. Thus far it 


has been determined what must be in the eset, all translational grounds. It 
must now be determined what else is allowed in the eset. 


Consider the zero internal strain energy requirement: 


a 


Taking the first row of equation (4.20b) following matrix multiplication 
yields: 
[KE ]{¥.} = KEI} (4.25a) 


Repartitioning this equation into translational and rotational dofs: 


ee tea ey th who! [Keols T + (4.25b) 


ie [eeu enya Keates [Ko] it tote 
Case III - if eset includes every translational dof and nog rotational dofs, 


equation (4.25b) becomes: 
Pa” oti Notte) 


43 


After matrix multiplication, equation (4.26) requires that: 


[Ke fer {te}s = {0} (4.27) 
which means that for every row, the sum of the translational dof columns for 
an unrestrained structure equals zero. This condition is always true. 
Therefore, this shows that all translational dofs are allowed to be included in 
the eset. 

Case IV - if eset includes every translational dof and one rotational dof, 


equation (4.25b) becomes: 


Hee esha fh? et ae 428) 


[Keler! [Kelee |ttele) = [9 [Keo}en }l {9} 
After matrix multiplication and the fact that Ke ee = {0}, equation (4.28) 
requires that: 
[Kel ellete = {0} (4.29) 
Since there is only one rotation in the eset, the dimensions of [Ke = are the 
# of eset translations by 1. Consequently, in order for equation (4.29) to be 


true, [Keele must always equal zero. By the definition that this case includes a 


rotational dof value in the eset, [Kole +[0]. Therefore, one rotational dof 


may not be included in the eset. 
Case V - if eset includes every translational dof and greater than one 


rotational dof, the same requirement (equation 4.29) as stated in Case IV 


exists. However, for this case, [Kole does not have to identically equal zero, 


44 


but the sum of the columns for every row of [Ke = must equal zero. In 


general, this sum is not zero. Therefore, it is concluded that no rotational dofs 
can be included in the eset. 
Case VI - if eset includes a subset of all translational dofs, equation 
(4.25b) becomes: 
ten hick E {eh [Kel Key (4.30) 
O Oll{ld. 0 0 || {0} 
After matrix multiplication and rearrangement, equation (4.30) requires that: 
[Ke {1-}, +[Ke]_ {le}, = {0} (4.31) 
Equation (4.31) is the same as if summing all translational columns of an 
unrestrained structure for every dof in the eset. From Case III (equation 4.27), 
it is seen that this sum is zero and equation (4.31) is satisfied. Therefore, this 
shows that any subset of all the translational dofs may be included in the eset. 
As a result of Cases I-VI, the following requirements on the choices of 
eset dofs for static displacement synthesis are summarized: 
(1) eset must include all translational grounded dofs 
(2) structure may be grounded anywhere 
(3) any subset of, or all translational dofs may be included in the eset 
(4) no rotational dofs are allowed in the eset 
For the work to be presented in this thesis, the requirements necessary for the 


reduced force vector to be determined from the reduced mass matrix are met. 
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Therefore, the task is now to derive a methodology for extracting [M] from 


FRF data. 


3: Formulation of Reduced Mass Matrix [M’] 
Since [H..(0)] =[K]"', it follows that: 
[H..(Q)] = [K-Q?Mr" (4.32) 
Taking two derivatives of equation (4.32) with respect to Q yields: 


d?H,,.(Q) _ 
da? 





=(K- oPNAT 120%) = ec m= +(Ik- Q?M]"[2M] + aa “= (2084) IK -0°N4)" 
(4.33) 


Applying the static condition of Q=0 to equation (4.33), and solving for [M1]: 


d°H,.(0) = 


40? (K] (4.34a) 


~ eee 
M]=—[K 
[M] al ae 
or rewritten in terms of the synthesis process: 


-1 d?H?:(0) 


same 
IM']==[H.@] 4 


[H.(0)] (4.34b) 
From equations (4.34a&b), the unknown second derivative of the eset FRF at 
zero frequency is present. The calculation of this value will be handled using 


the following forward finite differencing scheme with error order O(AQ?”): 


— 2 ( AQ) (4.35) 


The accuracy of this calculation is very sensitive to the size of AQ. The proper 


choice of AQ is dependent on the relative magnitudes of the mass and 
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stiffness matrices. It was found for this work that AQ = 0.1 provided a very 


dH, (0) 


ce However, this value is not universal and the 


accurate estimation of 


proper selection of AQ should be based on individual structures. 
From the previously discussed requirements on the conversion of [M] 


to {F}, it can be seen that static displacement synthesis is not arbitrary in the 
selection of which dofs may be included in the eset. However, this restriction 
does not apply to which dofs can be chosen to be modified, but only to 
whether static displacement information can be obtained about these dofs. For 
example, modifications to rotational dofs may be made to a structure and the 
synthesized FRF obtained. As a result of this, rotational dofs exist in the eset, 
which is not allowed. The solution to this dilemma is to re-synthesize the 
structure with no changes and only translational dofs in the iset. What this 
does is remove the rotational dofs from the eset and condense the FRFs prior 
to using them to obtain the reduced stiffness matrix and force vector. 
Therefore, the static displacement method is really not restrictive in its 


application, but rather in what information it can provide. 


D. STATIC DISPLACEMENT SYNTHESIS COMPUTER CODE 
The goal in the formulation of the static displacement synthesis 
computer programs is to first perform the frequency domain synthesis for an 


indirect modification operation to a given structure. The program allows for 
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the same type of modifications to the structure that were presented in the 
frequency domain synthesis chapter of this thesis. The program then uses the 
synthesized FRFs and returns the static displacements for all dofs designated 
in the eset. 

The reason why the frequency domain synthesis is performed here 
seperately vice using the FRFs obtatined from doing a frequency domain 
synthesis problem is two fold. First of all, there is no damping effect for a 
static problem. Therefore, the synthesized FRFs must be generated with zero 
damping. Secondly, if we want to determine the static displacement for a base 
excitation problem, the base excitation form of frequency domain synthesis is 
not what is desired. The correct form of the frequency domain synthesis is the 
force excitation synthesis equation, with the base excitation interconnection 
elements considered as springs to ground. Furthermore, an entire spectrum 
of FRF information is not necessary for static displacement calculations. 
Primarily only the static condition of Q=0 is needed. However, in this 


implementation, the synthesized FRFs at the first four frequencies are needed 


2+ 1* 
for the finite differencing calculation of a 
Appendix B shows the computer codes used to perform the static 
displacement synthesis, and perform a comparative analysis of the synthesis 


versus classical Guyan reduction methods. The programs presented in 


Appendix B can be seperated into two categories. The first category is a 
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program which uses the classical Guyan reduction techniques to calculate 
static displacement. The second category is the calculation of the synthesized 
static displacements where the different aspects of the synthesis technique 
have been modularized into MATLAB functions. The following is a diagram 


of the flowpath of the modularized static displacement synthesis programs: 






Modifications Form Impedance 






statchange.m statdeltam 


Figure 4.1 Flowpath of the Modularized Static Displacement Synthesis 


Programs 


It can be seen that the flowpaths for static displacement and frequency domain 
synthesis are very similar. This is done intentionally so that all synthesis 
techniques can be driven by the same modifications to a structure. A more 


detailed description of each function is contained within Appendix B. 


E. ILLUSTRATIVE EXAMPLES 


The purpose of this section is to illustrate the use of the static 


displacement synthesis technique and validate its accuracy. The mass-plate 


49 


structural system used in the frequency domain synthesis example, along 
with the changes that were made to the structure, is also used here to 
illustrate static displacement synthesis. The following tables are comparisons 
of the reduced matrices, vectors, and static displacements calculated using 


frequency synthesis and Guyan reduction. 
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1.0e+04 * 


1 Los 
-0.0000 
-0.0000 
0.2276 
-0.3794 
0.0000 


-0.0000 
1.0000 
-0.0000 
-0.0000 
-0.0000 
-0.0000 


[K*]=[H2,(0)]” 


-0.0000 
-0.0000 
1.0000 
-0.0000 
-0.0000 
0.0000 


0.2276 
-0.0000 
-0.0000 
1.3415 
-0.5691 
0.0000 


-0.3794 
-0.0000 
-0.0000 
-0.5691 

1.4985 
-0.5500 


0.0000 
0.0000 
0.0000 
0.0000 
-0.5500 
0.5500 


Table 4.1a Static Displacement Synthesis Reduced Stiffness Matrix 


1.0e+04 * 


1.1518 
0.0000 
-0.0000 
0.2276 
-0.3794 


0.0000 
1.0000 
-0.0000 
0.0000 
-0.0000 


[K] =[T]'(K][T] 


-0.0000 
-0.0000 
1.0000 
0.0000 
-0.0000 


0.2276 
0.0000 
-0.0000 
1.3415 
-0.5691 


-0.3794 
-0.0000 
0.0000 
-0.5691 
1.4985 
-0.5500 


Table 4.1b Guyan Reduction Reduced Stiffness Matrix 


Oo GO © 


-0.5500 
0.5500 


d?H°.(0)  —H?,(Q,)+ 4H?,(Q,) - 5H2,(Q,) + 2H", (0) 


dQ? (AQ)’ 

1.0e-07 * 

0.0671 0.0294 0.0294 0.0271 0.1094 0.1470 
0.0294 0.0588 0.0147 0.0294 0.0593 0.0593 
0.0294 0.0147 0.0588 0.0294 0.0593 0.0593 
0.0271 0.0294 0.0294 0.0774 0.1414 0.1978 
0.1094 0.0593 0.0593 0.1414 0.3567 0.5049 
0.1470 0.0593 0.0593 0.1978 0.5049 0.8242 


Table 4.2a Static Displacement Synthesis 2nd Derivative Matrix 


d?H_.(0) 


——s— = [K]I2M][K 
qo” [K][2M]IK] 
1.0¢e-07 * 
0.0671 0.0294 0.0294 0.0271 0.1094 0.1470 
0.0294 0.0588 0.0147 0.0294 0.0593 0.0593 
0.0294 0.0147 0.0588 0.0294 0.0593 0.0593 
0.0271 0.0294 0.0294 0.0774 0.1414 0.1978 
0.1094 0.0593 0.0593 0.1414 0.3568 0.5049 
0.1470 0.0593 0.0593 0.1978 0.5049 0.8242 


Table 4.2b Guyan Reduction 2nd Derivative Matrix 


2 


0.1928 
0.0903 
0.0903 
-0.0492 
0.0422 
-0.0000 


[M° 1==[H. (o)] 


0.0903 
0.2940 
0.0735 
0.0620 
0.1417 
-0.0000 


0.0903 
0.0735 
0.2940 
0.0620 
0.1417 
-0.0000 


1 ee 


-0.0492 
0.0620 
0.0620 
0.1537 
-0.0095 
-0.0000 


—=—[H:.(0)] 


0.0422 
0.14177 
0.1417 
-0.0095 
0.4215 
-0.0000 


-0.0000 
-0.0000 
-0.0000 
-0.0000 
-0.0000 
0.2588 


Table 4.3a Static Displacement Synthesis Reduced Mass Matrix 


0.1928 
0.0903 
0.0903 
-0.0492 
0.0422 


0.0903 
0.2940 
0.0735 
0.0620 
0.1417 


[M] = [T][M][T] 


0.0903 
0.0735 
0.2940 
0.0620 
0.1417 
-0 
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-0.0492 

0.0620 

0.0620 

0.1537 

-0.0095 
0 


0.0422 

0.1417 

0.1417 

-0.0095 

0.4215 
0 


Table 4.3b Guyan Reduction Reduced Mass Matrix 


o 2 © © 


0 
0.2588 


Synthesis: {F*}=[M‘]{1, }g 


-141.6000 
-255.6080 
-255.6080 
-84.5960 
-285.0197 
-99.9984 


Guyan: 


(F}=[T]' {F} 


-141.6012 
-255.6102 
-255.6102 
-84.5966 
-285.0226 
-100.0000 


Table 4.4 Reduced Force Vectors 


Synthesis: {x,}) =[K*]{F‘} 


-0.0296 
-0.0256 
-0.0256 
-0.0316 
-0.0714 
-0.0895 


Guyan: {x,} 


ZK 
-0.0296 
-0.0256 
-0.0256 
-0.0316 
-0.0714 
-0.0895 


Table 4.5 Static Displacements (in.) 
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From Tables 4.1 - 4.5, it can be seen that the static displacement synthesis 
calculations match those of the Guyan reduction at every step of the process 


on its way to calculate static displacement. 
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V. TIME DOMAIN STRUCTURAL SYNTHESIS 


The following background information on time domain structural 
synthesis is taken from reference [8]. This portion of the thesis will outline 
the methodologies and formulations of the time synthesis technique and the 
classical techniques used to check its accuracy. There will also be a time 
comparison study done in order to quantify the computational efficiency of 
the synthesis technique. 

As mentioned previously, time domain structural synthesis is a 
methodology whereby changes can be made to a given structure, and its new 
transient response determined as the solution of a nonstandard 
nonhomogeneous Volterra integrodifferential equation (VIDE) of the second 
kind. This integral equation is the result of the combination of the impulse 
response functions (IRFs) of the baseline structure, the reaction forces 
generated from the modification to the structure, and the total dynamic 
response of the original system in terms of the convolution integral. This 
greatly differs from the classical method of evaluating a structural change. In 
the classical method, the basic elements which define the structure are 
reconstructed, and then a modal superposition, convolution integral or direct 
integration technique is used on the new structure to solve for the response. 

The motivation for the development of the time domain structural 


synthesis technique stems from the advantages and flexibility that is enjoyed 
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through the use of the frequency domain structural synthesis technique. The 
time domain synthesis can be divided into the same two classes of 
modification and coupling, direct or indirect. Some of the advantages and 
flexibility shared by the two methods are: 1) the formation of an exact 
solution, with the ability to handle damping modifications; 2) arbitrary model 
reduction in the sense that only information about those dofs which are 
needed is retained; 3) the solution phase of the finite element analysis is only 
done once, and the new synthesis model elemental or system matrices are 
never formed. 

Just as with the frequency domain synthesis, this thesis will only focus 
on indirect modifications to a given structure. It will also allow for the 
inclusion of a spring stiffener or a viscous damper between two points, a 
lumped mass on the structure, and implementing a base excitation to the 
structure. A generalized computer program will be presented, which will 
allow for any of the above mentioned changes to be accomplished on some 
baseline structure. The program then returns the transient response for all 
dofs where changes have occurred and any other dofs that may be of interest 


to the user. 


A. GENERALIZED TRANSIENT RESPONSE 
The following equations are for determining the total dynamic 


response of a system which experiences some sort of excitation. These 
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formulations are taken from references [1 & 2], and are used to verify the 


accuracy of the synthesis after a modification has been made. 


il; Convolution Integral Method 
The convolution integral equation for determining the total dynamic 


response of a system is as follows: 


{x(t)}= {x(t}, + [ne —t)|{F(r)}de (5.1) 


where: 


{x(t)},, - homogeneous solution which contains the constants of 
integration 


h(t) - matrix of impulse response functions (IRFs) 
F(t) - excitation force 
t - dummy time variable 

The IRF matrix is determined modally using the equation: 


—_ @t 
Ch(t)] =() —~e oy r sin Oat [oy (5.2) 
d 


r 


where [®] = the assembled matrix of mass normalized mode shapes. Any 


element of the IRF matrix may be represented as: 


nmodes 1 —¢ @ t 
h..(t) = ‘¢} ——e Tf !l sin@a_t (5:5) 
g(t) re : 
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This is a relatively efficient means of calculating a structure’s dynamic 
response. However, this method is based on the principal of superposition 
which is valid for linear systems only. Also the modal method of obtaining 
the IRFs does not readily handle non-proportional damping. Therefore, if a 
damping change to the structure is made, or if a non-linear interconnection 
element is used (time domain synthesis with non-linear elements is possible, 
but will not be explored in this study), the response cannot be determined 


using the convolution integral method. 


ae Direct Integration Method 

The direct integration method of calculating a system’s dynamic 
response is able to handle the shortcomings experienced by the convolution 
integral method. However, the price is paid in the computational time 
necessary for this solution. For the work illustrated in this thesis, MATLAB’s 
ODE45.m function is used for this calculation [Ref. 9]. In order to use this 
function, the 2nd order ordinary differential equations (ODEs) of motion: 

[M]ix} +[C}ixt +[K]ix} = {Fy} (3.2) 

must be rewritten as a system of Ist order ODEs. The following is the general 


matrix equation which converts any linear system’s equations of motion to a 


Ist order system. 
lie} : y va) (ellie 3 heal (5.4) 


where: 
{z,} = {x} {z,} = {x} 
{Z>} = {x} 12> } = {x} 


ODE45.m solves for the vector [{z,} {z,}]" which are the system’s response 


displacements and velocities using an adaptive time stepping procedure. 


B. TIME DOMAIN SYNTHESIS FORMULATION 
Consider again the exact same general ndof structure which was 
introduced in Chapter 3 / Section B, and to which structural modifications are 


to be made: 


a 


Figure 5.1 General NDOF Structural System 


Using the convolution integral, equation (5.1), and the previously 
defined sets, iset, cset, bset, zset, eset for partitioning, the total dynamic 
response of the system can be written as: 


[ha(t-)] 1 [h.(t—2)] | [hat-a)] FEC) ¢dt 5.5) 


a f t[[hu(t—2)] | [he (t—2)] } fhylt— 2} ]/ RO 


Xe) (Xp), Oe ele Th, ¢— t)|: [h,, (t-7)] | F,(2) 

Before modification, the force vector refers to externally applied forces 
to the structure. However, following modification, there exists coupling 
forces between the cset and bset dofs. By definition, the iset dofs do not 


experience these coupling forces. Therefore, the partitioned forces may be 


rewritten as: 


F =F* (5.6c) 
F. =F + F° = F™ —(AM.X; + AC.x? + AK.x®) (5.6b) 
F, = Fy“ + F, =F“ -[AC,(%; — y) + AK, (x; — y)] (5.6c) 


where: 

f* - the coupling forces due to structure modification. 

[AM], [AC], [AK] - the modification matrices 

x", x", xX, - generalized (synthesized) responses after modification. 
Substituting equations (5.6) into equation (5.5), and recognizing that the 


presynthesized responses {x, x, x,}" are those due to the external forces only, 
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the following result is obtained: 


t 
fe a ffecsae esol Shed lS he FS tales 
x} |x} J olffit—2]" fh, DT vay ble ~Y 
(5.7) 
Equation (5.7) is the nonstandard nonhomogeneous VIDE of the second kind. 
In order to obtain the solution for this equation, the integral is discretized 
using a trapezoidal quadrature rule. Also by assuming that there are no 


external excitations to the structure, {Fo} = 0, equation (5.7) becomes: 


x) [TA] | [Aw] 
xf fie tall Pe} fee fer }) ea 
emi elet (5.8b) 


The elemental form of the quadrature matrix is 


(hy) 0 0 0 0 
Des 21.0 

A At 
A 2 (s) S(hj), 9 9 0 A 
[As] “n,), at(h,), (rs), 9 9 | 

_ ; 
A 
= (hs), Auh,)  At(h,) | --- (ny), 


where the subscript n represents the number of time steps, and the coupling 
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forces equal: 


aes shat (5.10a) 
i 7 - ection , (5.10b) 


tre fo[ oO a ¥ (5.100) 


e 


Equation (5.8b) can be rearranged and placed in the general form of 
[A] {x*} = {b} where {x"} can be solved for directly using various methods. 
However, due to the possible large size of [A], it is advantageous to leave 
equation (5.8b) as is, and solve for {x’} by iteration. The iteration procedure is 
as follows: 

1) assume the force vector {F,} 

2) calculate the response vector {x,}° 

3) use {x,}° to calculate {x,}", {x,}’,and a new {F,} 

4) use the new {F,} and calculate a new {x,}* 


5) repeat steps 3 & 4 until {x.} new - {X-} og S Convergence tolerance 


Although the time domain structural synthesis is exact in its 
formulation, approximations have now been introduced through the 
solution techniques of the trapezoidal quadrature, and the iteration method. 


Not only is the convergence important, but the stability of this solution is also 


of concern. The stability and convergence is dependent on the norm of A,, 
and the magnitude of the modifications AM, AC, AK. The norm is defined as 
the largest singular value of the matrix and in this case is driven by the value 
of At. The smaller the value of At, the better the stability and faster 
convergence is reached. The larger the value of the modifications, the greater 
the possibility for instability in the solution. Therefore, due to the solution 
technique, there are now limitations on the modification values. Inherent to 
and dependent on the computer system used, there will be limitations on the 
amount of time which can be covered by this analysis due to memory 


capabilities. 


Sc TIME DOMAIN SYNTHESIS COMPUTER CODE 

The goal in the formulation of the time domain synthesis computer 
programs is to perform the synthesis for an indirect modification operation to 
a given structure, allowing the same type of modifications presented 
previously in the frequency domain synthesis chapter of the thesis. The 
program then returns the transient response for all dofs where changes have 
occurred and any other dofs that may be of interest to the user. Appendix C 
shows the computer codes used to perform the time domain synthesis, and 
perform a comparative analysis of the synthesis versus classical methods. 

The programs presented in Appendix C can be segregated into two 


categories. The first is the calculation of the synthesized transient responses 
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due to a base excitation on the structure, and their comparison to the 
classically calculated responses. One of the programs in this category uses the 
convolution integral method in order to check the synthesis, while the other 
uses direct integration. The second category is the calculation of the 
synthesized transient responses where the different aspects of the synthesis 
technique have been modularized into MATLAB functions. The following is 
a diagram of the flowpath of the modularized time domain synthesis 


programs: 


Load IRFs 


change.m 


Load Structure Calculate IRFs Modifications 


change.m modal.m change.m 
impmodal.m 





Input Base 
Excitation 


fBlastForcing.m 





Build Quadrature Perform Synthesis 
Matrix 





buildA.m timesynth.m 


Figure 5.2 Flowpath of the Modularized Time Synthesis Programs 


A more detailed description of each function is contained within Appendix C. 
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D. ILLUSTRATIVE EXAMPLES 

The purpose of this section is to illustrate the use of the time domain 
synthesis technique and validate its accuracy. A computational time 
comparison will be done between the synthesis and classical methods in order 
to assess the computational efficiency of this technique. The three structural 
systems previously described, mass-spring-damper, beam and mass-plate, are 


used to accomplish this goal. 


1. Mass-Spring-Damper System 


Consider the following mass-spring-damper system: 


m, ey k = 1000 lb/in 
\ eee m = 0.259 lb-s?/in 
kX, ! Ak,, = 250 Ib/in 
—lAMllAAIF Ak, = 250 Ib/in 
ji- 2.9 m1 Am, = .02 Ib-s’/in 
CMe: ky Ac, = 2 lb-s/in 
y(t) ) ee r y(t) = Y,(e"" -e™) 
Se a 0 


Figure 5.3 Mass-Spring-Damper System Experiencing Base Excitation 


The solid elements represent the original structure and the dashed elements 
represent modifications. Although not shown in the figure, the original 
structure is proportionally damped using [C] = a[K]. The modifications shown 


are arbitrarily selected. 
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a. Spring Modification Only 

The following figure is the resultant synthesized and classically 
determined transient response x,(t) which represents the response at dof 2 
due to the base excitation. For this example neither the mass nor the damper 


modification was included. 


Transient Time Response at dof 2 







displacement (in) 


_ 0.1 
time (sec) 


Figure 5.4 Transient Response Mass-Spring-Damper System 


From the graph produced, it can been seen that the responses are identical, 
proving that the synthesis technique is exact, and that the computer coding 
for the arbitrary changes is correct. These results were produced using the 


tsynconv.m program. 


68 


b. Spring & Mass Modifications 

The previous example is now repeated but with the mass change 
at dof 2, Am, = .02 Ib-s*/in included. The following figure shows the results of 
this additional modification. 


Transient Time Response at dof 2 





displacement (in) 


0.1 


0.3 ! : 
0 0.05 _ 0.1 0.15 0.2 
time (sec) 


Figure 5.5 Transient Response Mass-Spring-Damper System w/Spring-Mass 


Modification 


From the graph produced, it can been seen that the mass modification had no 
effect on the accuracy of the solution. However, the addition of the mass 
modification did have a large impact on the computational time required for 


the synthesis method. This will be discussed in more detail later. 
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c Spring & Damper Modifications 
The spring modification only example is now repeated but with 
a damper change at the base Ac, = 2 lb-s/in included. The following figure 


shows the results of this additional modification. 


Transient Time Response at dof 2 
0.3 


0.25 


= 
ND 


© 
= 
mn 


2 
a 


0.05F- 


displacement (in) 


© 


-0.05 


-0.1 





-0.15 : 
0 0.05 


_ 0.1 
time (sec) 


Figure 5.6 Transient Response Mass-Spring-Damper System w/Spring- 


- Damper Modification 


From the graph produced, it can been seen that the responses are identical 
with respect to the resolution of the plot graphics. Also, to further validate 
the correctness of the computer algorithms, the effect of the added damper 
can be seen when comparing Fig. 5.4 & Fig. 5.6. Since the non-proportional 


damper change was made to the structure, the direct integration method was 
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used to obtain the classical solution. Although there was no real effect in the 
accuracy of the solution, the effect in computational time is profound for the 


classical method. These results were produced using the tsynode45.m 


program. 


d. Spring , Mass, & Damper Modifications 
The example is now repeated with all modifications included. 


The following figure shows the results of these modification. 


Transient Time Response at dof 2 


0.3 


0.25 


0.2 


0.15 


0.1 


0.05 


displacement (in) 


-0.05 


-0.1 





“0.15, 0.1 
time (sec) 


Figure 5.7 Transient Response Mass-Spring-Damper System w/Spring-Mass- 


Damper Modification 
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Just as in all the other cases, it can been seen that the responses are extremely 


close to identical. These results were also produced using the tsynode45.m 


program. 


Ds Beam System 


Now consider the following beam system: 


Ak, = 1000 Ib/in 


fREEE EEE FD Acyso = 2 Ib-sfin 
! y(t) = Y,(e% - €°) 


? 
(ee 
. k, cy : Yo= 
’ 


Figure 5.8 Beam System Experiencing Base Excitation 


The beam parameters and finite element modeling of the above structure is 
identical to that used in the frequency domain synthesis chapter of the thesis. 


The modifications shown are arbitrarily selected. 


a. Base Spring Modification Only 

The following figure is the resultant synthesized and classically 
determined transient response x,,(t) which represents the response at dof 11 
due to the base excitations. As can be seen from Fig. 5.7, this is not a dof where 


change has occurred or where the base excitation is attached, but just a dof 


72 


where there is 


interest in the response. For this example the damper 


modification Ac,,, was not included. 


displacement (in) 


Transient Time Response at dof 11 
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Figure 5.9 Transient Response Beam System 


From the graph produced, it can been seen that the responses are not 


identical, but follow each other with relatively little error. The difference 


between the two methods can be attributed to the numerical solution of the 


synthesis method. As the time step for analysis is decreased, the numerical 


solution approaches the exact solution. However, the use of smaller time 


steps is restricted by computer memory limitations. 


is 


b. Damper Modification 
The previous example is now repeated but with the damper 
change at dof 13, Ac,,, = 2 Ib-s/in included. The following figure shows the 


results of this additional modification. 


Transient Time Response at dof 11 
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Figure 5.10 Transient Response Beam System w/Damper Modification 
From the graph produced, it can been seen that the responses are not 
identical, but follow each other with relatively little error. The same line of 


reasoning as in the undamped case is applicable here also. A smaller time step 


yields more accurate data but, memory limitations inhibits its usefullness. 
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2. Mass-Plate System 


Now consider the following mass-plate system: 







Ak, = 1e4 Ib/in 
Ac, = 5 lb-s/in 
Ak, = 500 Ib/in 


Ac, = 5 lb-s/in 


y(t) = Y,(e™™ - 7) 
Y,=1 


Figure 5.11 Plate-Mass System Experiencing Base Excitation 


The plate parameters and finite element modeling of the above structure is 
identical to that used in the frequency domain synthesis chapter of the thesis. 


The modifications shown are arbitrarily selected, but are in accordance with 


the changes which will be made during an optimization problem. 


a. Spring Modifications Only 


The following figure is the resultant synthesized and classically 


determined transient response x,).(t) which represents the response at dof 109, 
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the computer, due to the base excitations at the plate’s corners. For this 


example none of the damper modifications are included. 


Transient Time Response at dof 109 
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Figure 5.12 Transient Response Mass-Plate System 


From the graph produced, it can been seen that the responses are not identical 
but follow each other with relatively little error. The error difference between 
the two methods can be attributed to the numerical solution technique, and 
the step size of the time vector. As the step size decreases, the synthesis 
solution approaches the exact solution. However, memory limitations on the 
computer system used for this research, prohibit meaningful use of a smaller 


time step. 
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b. Spring & Damper Modifications 

The spring modifications only example is now repeated but with 
the prescribed damper changes at the base and between the computer and the 
plate. The following figure shows the results of these additional 


modifications. 


Synthesized Transient Time Response at dof 109 
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Figure 5.13 Transient Response Mass-Plate System w/Spring-Damper 


Modifications 


Due to the size of this model and the damper modifications being made, it 
was not possible to compare the synthesis solution with a classical solution. 
The computational time required for a classical solution would be great. 


However, from comparing Figs. 5.11 & 5.12, the slight effect of the additional 
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damping can be seen in the peak amplitude of the response, giving some 
credence to the solution. Also, the other previously run damper modification 


examples show that this technique does in fact work. 


4, Computational Time Comparisons 
The following table is a compilation of all the computation times that 


it took for the previously presented examples to complete its process. 


Synthesis Classical Methods 


Convolution | Integration At 
| cy | es) | ee) | ea 
Mass- No Damper 15162 1.28 O01 
spring | Nomass | | | 
Damper Zale o/Foe 001 
No Mass | 


Mass 194.05 001 
No Damper a 

Mass 197.89 I39=22 001 
Damper | 

: 0003 


5687.88 0003 


Damper 183.75 


Beam No Damper 107.38 


Mass-Plate | No Damper 27.39 0008 
Damper Ig9 t > 00 0008 


Table 5.1 Synthesis vs. Classical Computational Times 
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From the table, it can be observed that for smaller structures, the classical 
method seems more efficient. This is due to the fact that the advantage of the 
iteration method is lost when dealing with small structures. The 
computational time for the iteration solution, or in other words the synthesis 
technique, is driven by the magnitude of the time step and the number of 
time steps. On the other hand, the classical methods’ solution times are 
driven by the model size. Therefore, depending on the time step size, a very 
small and very simple model may take a very long time to solve using time 
domain synthesis. 

Another interesting aspect of the comparison is how sensitive the 
synthesis technique is to mass changes. Since the mass change does not add 
additional dofs to the model, the classical method’s computational time is not 
increased by this modification. On the other hand, a mass change for the 
synthesis method constitutes the use of the calculated second derivative of 
the synthesized response. This derivative is accomplished using a finite 
differencing scheme and inherently adds the possibility of additional 
instability to the solution. 

As previously alluded to, the finite differencing could be the reason 
why the synthesis method is greatly affected by mass changes. However, 
damper changes, which also prompt the use of finite differencing, did not 
substantially increase computational time for the synthesis method. This is in 


contradiction to the thought that the increased instability is caused by the 
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finite differencing scheme. Another factor which plays into the convergence 
and stability of the solution is the magnitude of the modification. Relative to 
each other and this problem, the mass change may be of greater magnitude 
than the damper change and therefore costs more in computational time. 
The largest and most noticeable disparity is in the time required to 
classically calculate the dynamic responses following a damper modification. 
The direct integration method is very inefficient and becomes basically 
unusable for very large complex models, or in the use of an optimization 


routine. 
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VI. OPTIMIZATION 


The following background information and the motivation for this 
chapter can be attributed in part to reference [7]. Inherent to the operation of 
military vehicles, whether they operate on land, sea or in the air, is the 
presence of a shock and vibration environment. Computer and electronic 
equipment located on these vehicles can malfunction or be damaged as a 
result of the severe effects of this environment. Therefore, the need arises for 
a properly designed isolation system, which will absorb the shock and 
vibrations and therefore protect the equipment. This is an especially 
important concept in today’s military with the relatively recent practice of 
using Commercial-Off-The-Shelf (COTS) equipment vice the traditionally 
expensive shock and vibration hardened military equipment. 

In general, it is relatively straightforward to design an isolation system 
to protect against shock only or vibration only. However, these two designs 
are in direct conflict with each other [Ref. 1]. Therefore, in the presence of 
both types of excitations and various other constraints, the optimal design of 
an isolation system becomes more complicated. The need then arises for 
some sort of optimization scheme which will minimize the response of the 
component to be protected, while simultaneously satisfying the constraints 


which have been placed on the system. 
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The complexity of this problem does not necessarily end here. If the 
structure to be isolated and the isolation system to be designed is large and 
complex, the analysis process to calculate the desired responses can be very 
computationally demanding. Furthermore, if there are multiple design 
variables which may be altered in order to achieve the desired result, the 
search patterns for the optimal solution are more complex, which leads to 
more iterations of the optimization process. The combination of these two 
factors can lead to the conclusion that the optimal design of a large complex 
isolation system is impractical. 

With the previously presented arguments in mind, the use of 
computationally efficient synthesis techniques for the calculation of the 
system’s response is the obvious answer. Since the structure is large and 
complex, analysis time could be relatively long. The larger the system, the 
greater the number of possible design variables. This means that the search 
for the optimal solution is even more complex, and will require more 
iterations. This portion of the thesis will demonstrate the use of the 
previously presented structural synthesis techniques in the optimal design of 
a shock and vibration isolation system. A general computer algorithm will be 
formulated and used in order to illustrate the use of structural synthesis in 


optimization design. 
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A. GENERAL OPTIMIZATION PROBLEM FORMULATION 

In an optimization problem, the quantity to be maximized or 
minimized is termed the objective function. The quantities which may be 
altered in order to find this optimal value are termed the design variables. 
The constraint functions are also a function of the design variables, and three 
different types of constraints may exist in an optimization problem. The first 
type is an inequality constraint, where the value calculated from the 
constraint function is less than or equal to some prescribed limit. The second 
type is an equality constraint. Here the value calculated from the constraint 
function must be equal to some prescribed value. The third constraint type is 
side constraints. This is when upper and lower limits are placed on the design 
variables. In a constrained optimization problem, the objective function is 
maximized or minimized, while ensuring that the user defined values of the 
constraint functions are not violated. 

For the design of an isolation system, the minimization of the 
equipment’s response over some frequency range, due to vibration, could be 
the main concern. At the same time, the entire structure’s response due to a 
shock input, and the static displacement or ‘sag’ of the system could be 
constrained. The physical parameters which can be modified to achieve this 
goal could be the stiffness and damping values of the isolators. The design of 


this isolation system can be converted to the following optimization problem: 
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Minimize (Objective Function): 
Computer response due to random vibration 
Subject to (Constraints): 
Maximum static displacements < limit 
Maximum dynamic isolator stretch due to shock ¢ limit 
Maximum computer acceleration due to shock < limit 
Min and max changes in the isolator stiffness values < limit 
Min and max changes in the isolator damping values < limit 
Design Variables: 
Changes in the isolator stiffness values 


Changes in the isolator damping values 


1. Calculation of Objective Function 

Since the objective function is based on the response due to random 
vibrations, the input base excitation is in the form of a power spectral density 
(PSD) function. The response is therefore calculated as the mean square value 
in terms of the system response function (FRF), and the PSD of the input. The 


equation for the mean square value of the response is given as [Ref. 10]: 


se 2 
x’ =| lA(f)| S(fdf (6.1) 
From equation (6.1), it can be seen that the most accurate means of 


minimizing this equation would be if the input PSD, S(f), is known. Since the 
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goal is to formulate a general optimization program, it is not inconceivable 
and is appropriate to assume S(f) to be constant over some frequency, f,<f < 
f,. Therefore, minimizing the mean square response over the appropriate 


frequency range is equivalent to minimizing: 
fr 2 
F= | "(PY df (6.2a) 


It can be seen that for every iteration of the optimization process, the 
calculations involved in equation (6.2a) must be performed. If the FRFs of the 
baseline structure is provided or calculated prior to the iterations begin, the 
frequency domain structural synthesis technique may be used to calculate all 


subsequent FRFs and the objective function is properly written as: 


aca 
F= | "f(f) of (6.2b) 
Zs Calculation of Constraints 
The first constraints mentioned are the maximum static displacements 
of the structure. Since these values must also be recalculated every iteration, 
an efficient means of doing so needed. Various reasons make this problem a 
prime candidate for the use of static displacement synthesis: 

e the static displacement at select locations vice the entire structure is 
desired, therefore calling for some type of reduction or 
condensation technique. 

e the elemental and system matrices of the structure may not be 


available. 
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e the modal solution of the eigenvalue problem has already been 
accomplished. 

e the isolator stiffness, which are design variables, act in vertical 
translation direction only (isolator damping has no effect on static 
problem). 

Another constraint is the maximum dynamic isolator stretch due to 
the shock input. In order to calculate this constraint, the transient response of 
the structure at the dofs where the isolators are connected must be obtained. 
The displacement of the base input as a function of time must also be known. 
Since the change in the isolator damper represents nonproportional 
damping, the direct integration method would be used to solve this problem. 
This calculation would be terrible, in reference to computing time for a single 
iteration of the optimization process. Despite other reasons, this alone is 
enough motivation for the implementation of the time domain synthesis 
technique for the calculation of the dynamic responses. The last constraint, 
the computer acceleration, is obtained by simply applying a finite differencing 
scheme to the computer’s displacement time history. All the above 
mentioned constraints are known as inequality constraints, and will satisfy 


the general normalized form g,(X) < 0 [Ref. 11]. 


86 


3. Design Variables 

There are limits on the amount in which the isolator stiffness and 
damping can be changed. For example, the lower bound for the change in an 
isolator stiffness should not be equal to the negative of the original value of 
the isolator. If it is, the optimization program might use this value and this 
constitutes the total removal of the isolator and changing the basic structure 
of the system. The upper bound for the change is limited by what is physically 
realistic. These type of constraints on the design variables are known as side 
constraints [Ref. 11]. 

The relative value of design variables to one another is very important 
when considering the efficiency and reliability of the numerical optimization 
process [Ref. 11]. The function space of the problem, which is defined by the 
design variables, should be as symmetric as possible. When looking at the 
design variables, isolator stiffness and isolator damping, it can be seen that 
effective changes to isolator stiffness are on the magnitude of thousands, 
while the isolator damper changes are on the magnitude of ones. This large 
disparity creates a very nonsymmetric function space. In order to alleviate 
these problems, the design variables are scaled so that their magnitudes are 
comparable. A symmetric function domain allows for the faster 
determination of the optimization by decreasing the number of iterations 


necessary to find the optimal solution[Ref. 11]. 
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B. OPTIMIZATION COMPUTER CODE 

The optimization tool used to solve this problem is the MATLAB 
optimization toolbox 1.0d. More specifically, the function constr.m, which is 
designed for optimization of a constrained, nonlinear, multivariable 
problem, is used. This function uses the Sequential Quadratic Programming 
(SQP) method in order to solve this problem [Ref. 12]. The SQP method is 
used in order to determine search directions for the design variables at every 
iteration [Refs. 11 & 12]. One of the limitations of this, and other traditional 
nonlinear optimization codes, is that they may only give local solutions. 
Also, when the problem is infeasible, constr.m attempts to minimize the 
maximum constraint value. What this can lead to is, the problem iterating to 
the user defined iteration limit, the search pattern extending outside the 
bounds of the design variables, and no optimum solution ever being found. 
Therefore, the user must have good understanding of the problem and the 
constraints imposed upon it. 

The goal in the formulation of the optimization computer programs is 
to provide a means of determining the optimal changes in isolator values, in 
order to minimize some dynamic response, while simultaneously satisfying 
all prescribed constraints. These programs first allow for the loading of a 
general finite element model of the structure. This baseline structure is given 
in terms of its FRFs and IRFs, or in terms of its system matrices from which 


the baseline FRFs and Ifs must be calculated. Modifications to the baseline 
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structure which will act as design variables are then designated. User defined 
initial values of design variables and constraint values are designated. During 
the optimization process, all structure responses will be calculated using 
structural synthesis techniques. With minor coding adjustments, the ability 
to interchange objective function and constraints exists. For example the 
optimization program could be the minimization of the computer’s 
acceleration due to shock, while the response to a random vibration input is a 
constraint. Appendix D shows the computer codes used to perform this 
process. The following figure is a diagram of the flowpath of the optimization 


program. 
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Form Change 
Matrices 


delta.m 


Calculate FRFs & IRFs 


Calculate Static 
Displacement 


statsynth.m 


yss 


Display 
Results 


Input Shock 
Excitation 


fBlastForcing.m 





Calculate 
Constraints 


isosub.m 





Figure 6.1 Flowpath of the Optimization Program 
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Cc. ILLUSTRATIVE EXAMPLE 

The purpose of this section is to illustrate the use of the structural 
synthesis techniques in the optimization of a shock and vibration isolation 
system. The figure below shows the mass-plate structural system which was 


used in the previous examples and is also used here. 


% yk fy Ty OT 


Figure 6.2 Plate-Mass System Experiencing Base Excitation 


i Problem Formulation 
The following is the formal problem statement of the optimization 


problem: 


2 
Minimize: F(Ak, ,Ak,,Ac,,Ac.)= [- [E(p)| df 


Z,_Stat 


51 = max z,_Stat - 
g, = Z._Stat -1<0 
max Z._ Stat 
eee k,_stretch 
ubject to: ee eee 
63 max k ,_stretch 
__k. stretch _ 
Be max k._stretch 
] 
g. Z,_acce 1<0 


max Z._accel 


~4000 < Ak, < 5000 
~4000 < Ak. < 5000 


O< Ac, <5 
O<Ac. <5 


where: 

H"(f) = frequency response function of the computer 

z, Stat = plate static displacement [in] 

zZ,_stat = relative static displacement between computer and plate [in] 

maxz, stat = maximum allowable plate static displacement [in] 

maxz,_stat = maximum allowable static displacement between 

computer and plate [in] 
stretch = isolator stretch between base and plate due to shock [in] 
k, stretch = isolator stretch between plate and computer due to shock 
[in] 

max k, stretch = maximum allowable isolator stretch between base and 
plate [in] 

max k,_stretch = max allowable isolator stretch between plate and 
computer [in] 

z,_accel = absolute acceleration of computer due to shock [in/s’] 

max Z,_ accel = max allowable absolute acceleration of computer [in/s’] 

Ak,, A k,, Ac,, Ac, = the design variables. 


a2 


The following is information that must be provided to the 
optimization program for its execution, and for the understanding of the 
user. The following table lists the initial, lower bound, and upper bound 


values for the design variables. 







initial value 
lower bound 
upper bound 


1d = 


Table 6.1 Optimization Inputs 


From the table it can be noticed that an initial value for the design variables is 
used. Since the optimization routine locates local minimum, the start 
position of the search is very important. The baseline values for the isolator 
elements which are to be modified are: 
initial k, = 10000 Ib/in 
initial k. = 5000 Ib/in 
initial c, =0 Ib-s/in 
initial c, =0 lb-s/in 
The constraint values for the inequality constraint are: 
maxz,, Stat = -0.08 in 
maxZ, stat = -0.03 in 
max k, stretch = 1.0 in 


max k,_stretch = 1.0 in 
max z,_accel = 30¢’s 
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Ze Optimization Results 
The previously described problem is executed using the isomain.m and 
isosub.m programs. During the execution of the optimization program, 


output information is provided through the MATLAB Command Window 


(Table 6.2). 

f-COUNT FUNCTION MAX{g} STEP Pr ur 

5 36012.1 1.40603 I 

10 15494 0.0364448 I 

1S 6969.56 -0.0262839 1 

De 6691.27 -0.0370827 0.25 

29 6372.67 -0.0682749 0.25 

34 6372.59 -0.0693421 1 

39 6370.5 -0.0631356 i mod Hess(2) 
44 6370.49 -0.0687934 1 mod Hess 

49 6370.49 -0.0689084 1 mod Hess(2) 
54 6370.48 -0.067107 I 

59 6370.47 -0.0527099 1 mod Hess 
64 6370.45 -0.0667607 I 

69 6370.44 -0.0687255 1 mod Hess(2) 
74 6368.61 -0.0685718 1 mod Hess(2) 
79 6368.61 -0.0688113 I mod Hess 
80 6368.61 -0.0667838 il mod Hess 


Table 6.2 Optimization 1.0d Generated Output 


Once the optimization program terminates successfully, post 


processing occurs and the following figures and summary information is 


provided. The first figure is a graph of the value of the design variables at 


each iteration. 
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Change in Variables vs. Iterations 
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Figure 6.3 Values of the Design Variables 


The next figure is a graph of the absolute value of the isolators at each 


iteration of the optimization process. 
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Isolator Stiffness (Ib/in) 
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Change in Isolation Constants vs. Iterations 
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Isolator Damping (Ib-s/in) 
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Figure 6.4 Absolute Values of the Isolator Elements 


The next figure is a graph of the FRF of the computer before the optimization 


begin, and at its termination. 
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Optimized Frequency Response Function [H*] 
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Figure 6.5 Computer FRF Before and After Optimization 


The effect of the changes to the structure, particularly the damping additions, 
can definitely be seen in this figure. The objective function is defined as the 
area under this curve after it is squared. Therefore, by lowering the peak 
values of the response, the objective function is ultimately lowered. The next 
figure is a plot of how the value for the objective function changed as a 


function of the optimization iterations. 
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Change in Objective Function vs. Iterations 
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Figure 6.6 Values of the Objective Function 


From Fig. 6.6, the dramatic effects of the optimization can really be seen. The 
initial value of the objective function was 36012.1 and after the isolator 
modifications the objective function equals 6368.31. The following 
information is a summary of the optimization process and is displayed at the 


MATLAB Command Window: 


Optimization Terminated Successfully 


Active Constraints: 
ans =3 


The recommended change in base isolator stiffness (Ib/in) is ---> 
final_del_kb = 2.3271e+03 Ib/in 
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The recommended change in computer isolator stiffness (Ib/in) is ---> 
final_del_kc = 4.0849e+03 lb/in 


The recommended change in base isolator damping (Ib/in/s) is ---> 
final_del_cb = 2.7678 lb-s/in 


The recommended change in computer isolator damping (lb/in/s) is ---> 
final_del_cc = 4.6266 lb-s/in 


Constraints: 
g, = -0.1814 g, = -0.6334 g3 = -0.5237 g, = -0.6981 g, = -0.0668 


elapsed_time = 9.5550e+03 secs = 2hrs 39min 


It can be seen from the presented information that it took the 
optimization program 80 iterations to recommend the presented changes to 
the structure in order to optimize its response against random vibrations, and 
meet the established constraint criteria. The most notable portion of this 
example is that it accomplished this task on a 109 dof structure, which 
underwent non-proportional damping changes at every iteration, in only 2 
hours and 39 minutes. In other words, this program performed 80 transient 
response analyses, calculated the static displacement 80 times, and calculated 
the frequency response function and the integral of that curve 80 times, for a 


109 dof structure. 
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VIl. CONCLUSIONS / RECOMMENDATIONS 


The main objective of this thesis was to investigate the use of 
structural synthesis techniques as rapid reanalysis procedures in the optimal 
design of large shock and vibration isolation systems. In order to use these 
synthesis techniques liberally, an investigation into the accuracy of the 
methods, and the computer algorithms formulated, had to be accomplished. 
It is important to mention at this point that the time and frequency structural 
synthesis methods used in this study, primarily covers base excitations. This 
is a new extension of the more familiar synthesis formulations. The static 
displacement synthesis, is however a completely new development. These 
formulations had to be general enough to handle different structures and to 
handle various types of modifications to the structure. The synthesis 
techniques were then coupled with an optimization routine, and together 


they were used in finding the optimal changes to a given isolation system. 


A. CONCLUSIONS 
From the various analyses conducted, the following conclusions can be 
drawn: 
1. The frequency domain structural synthesis technique is a very 
efficient and accurate method of calculating the new FRFs of a 


structure following indirect modifications. It is nonrestrictive and 
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arbitrary as to where and what type of modifications can be made to 
the structure. 

2. The static displacement structural synthesis is also a very efficient 
and accurate method of calculating the new static displacement of a 
structure following modifications. The accuracy of this method does 
depend on the use of a finite differencing scheme which in turn is 
dependent on the variable step size. It is also found that this method 
is restrictive in its use. Certain criteria, which are identified in 
Chapter IV of this work, must be met in order that this synthesis 
technique can be applied properly. For the isolation problems 
presented in this work, the required criteria were met. 

3. The time domain structural synthesis technique is an efficient and 
accurate method of calculating the new transient response of a 
structure following indirect modifications. It is nonrestrictive and 
arbitrary as to where and what type of modifications can be made to 
the structure. However, since numerical techniques are used to 
solve for the responses, some error of approximation is entered into 
the solution. As the time step is decreased, the numerical solution 
approaches the exact solution. However, current solution 
algorithms are restricted by memory requirements. 

4. In finding the optimal design of a relatively large shock and 


vibration isolation system, the use of the synthesis techniques 
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proved to be an invaluable asset in the optimization routine. With 
the use of the synthesis techniques, in 2 hours and 39 minutes, it 
was possible to perform 80 FRF, static displacement, and transient 
analyses on a 109 dof system, which was non-proportionally 
damped. To get an idea of how efficient this is, recall from the time 
domain synthesis analysis (Chapter V), that it took approximately 1 
and 1/2 hours to calculate one transient analysis of a 22 dof beam 


structure. 


B. RECOMMENDATIONS 
Throughout this study, some weaknesses of the techniques, and 
opportunities to enhance the techniques were touched upon. 
Recommendations to address these weaknesses, and take advantage of these 
opportunities include: 
e All modal analyses were done to include all modes of the structure. 
A modal truncation criteria should be established, and studies 
should be done to ensure that these criteria hold up for synthesis 
also. This way it would be possible to take advantage of time saving 
aspects which are inherent to modal superposition techniques. 
e Use actual FRF and IRF data from a real structure, and predict 


structural responses due to modifications. 
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Build CAD models to validate damping modifications to large 
structures. 

Investigate to obtain a general rule for choosing the optimum AQ to 
be used in finite difference scheme to calculate static response. 
Another option would be to investigate the use of the closed form 


aT tT 
2nd derivative of H,, to obtain a 

dQ 
Implement relaxation schemes for the iterative solution in time 
domain synthesis. Also a study should be done in order to 


determine what the proper time step size should be to ensure 


convergence and stability. 
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APPENDIX A. FREQUENCY DOMAIN SYNTHESIS COMPUTER 


CODES 


The following is a list and a brief description of the main MATLAB 
computer codes that were written in order to perform the frequency domain 
synthesis calculations. 

e fsynstr.m - performs the frequency domain synthesis on a structure which 
experiences an external force excitation directly on the structure, and 
compares this solution to the FRF calculated using classical methods. 

e fsynbase.m - performs the frequency domain synthesis on a structure 
which experiences a base excitation, and compares this solution to the FRF 
calculated using classical methods. 

e fsynmain.m - the main program which calls the modularized frequency 
synthesis programs and does the post processing of the results. 

e change.m - loads the baseline structure to be modified, allows for 
modifications to the baseline structure, and calls the function delta.m 

e delta.m - forms the modification matrices which will be used to form 
impedance matrices or determine coupling forces. 

e modal.m - solves for the structure’s natural frequencies, mode shapes and 
other modal matrices and vectors. 


e frfmodal.m - calculates the structure’s FRFs. 
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e frfsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized FRFs. 

e frf sift.m - selects the FRF that the user wants to perform post processing 
on. 


The full codes are contained on subsequent pages. 
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% This program is designed to calculate FRFs by inverting the impedance matrix. 

% Changing the mass, stiffness, and damping matrices, new FRFs will be calculated by 
% reformulating the matrices, and using the synthesis method in the frequency domain. 
%o 


clear 


% Formation of the elemental matrices 
% 
nel=4; 
=nel+]; 
dof=nel+1; 


ndk=1; 
knel(1)=0; 
for spring=1:ndk; 
kcon(spring)=100; 
knel(spring+1)=4; 
for ele=knel(spring)+1:knel(spring+1) 
kcon_ele(ele)=kcon(spring); 
end 
end 


ndm=1; 
mnel(1)=0; 
for mass=1:ndm; 
m(mass)=100/386.4; 
mnel(mass+1)=5; 
for ele=mnel(mass)+1:mnel(mass+1) 
mele(ele)=m(mass); 
end 
end 


% Formulation of global stiffness, damping and mass matrices from elemental matrices 


Oo 
Kglo=zeros(dof,dof); 
Cglo=zeros(dof,dof); 
Mglo=zeros(dof,dof); 
for n=1:nel 
kele=kcon_ele(n)*[ 1 -1 
rcl=n; : 
rc2=n+1; 
Kglo(rcl:rc2,rc1:rc2)=K glo(re 1:rce2,rc1:rc2) + kele; 
end 
Mglo=diag(mele); 
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alpha=0; beta=.0005; 
Cglo=alpha*M¢glo + beta*Kglo; 


ea ze 

if pc——| 
K=K glo(2:dof,2:dof); 
K((dof-1),:)=[]; 
K(:,(dof-1))=[]; 
M=Mglo(2:dof,2:dof); 
M((dof-1),:)=[]; 
M(:,(dof-1))=[]; 
C=Cglo(2:dof,2:dof); 
C((dof-1),:)=[]; 
C(:,(dof-1))=[]; 
ndof=dof-2; 

elseif bo==2 
K=K glo(2:dof,2:dof); 

=Mglo(2:dof,2:dof); 
=Cglo(2:dof,2:dof); 

ndof=dof-1; 

elseif bc==3 
K=K glo(2:dof,2:dof); 
K((dof-1),:)=(]; 
K(:,(dof-1))=[]; 
M=M¢lo(2:dof,2:dof); 
M((dof-1),:)=(]; 
M(:,(dof-1))=[]; 
C=Cglo(2:dof,2:dof); 
C((dof-1),:)=[]; 
C¢:,(dof-1))=]; 
ndof=dof-2; 

else 
K=Kglo; §©M=Mglo; C=Cglo; 
ndof=dof; 

end 


% This portion of the program is designed to accept changes to the original 

% structure and recalculate the FRF using classical method of reforming [M], 
% [K], and [C] matrices using the inverse of the impedence matrix. Also 

% recalculates the FRF using structural synthesis in the frequency 

% domain. 

Jo 


change_model=1; 
while change_model==1; 


K.c=K;. (Micc=Mi gs Gic—G- % initializes change matrices to 
% original matrices 


changek=2; changec=2; changem=2; ~% initializes logic variables 
% to ‘no’ status 
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chk=0; chc=0; chm=0; % initializes counter for # of changes 
% to matrices 


kdofl=[]; kdof2={]; % initalizes the matrices which keep 
cdofl=[};  cdof2=[]; % track of the dofs where changes occur 
mdof=[]; 


cset=0; % initializes cset matrix to zero to avoid null index error in 
% handling the spring to ground if no changes are made to the 
% structure. 


changek=input(‘Would you like to change your stiffness matrix? l=yes 2=no '); 
while changek== 
chk=chk+1; 9% counter for determining # changes to stiffness matrix. 
1k(chk)=input(‘input value of added spring (lbs/in) '); 
kdof1(chk)=input(‘input constrained dof where Ist end of added spring is applied '); 
kdof2(chk)=input(‘input dof where 2nd end of added spring is applied or 0 for ground 
") e 


? 


% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 


O 
if kdof1(chk)~=cset 
cset=[cset kdof1(chk)]; 
end 
if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 


% Classical method of reformulatting the [K] matrix. 


K_c(kdof1(chk),kdof1(chk))=K_c(kdof1(chk),kdof1(chk)) + 1k(chk); 
if kdof2(chk)~=0 
K_c(kdof1(chk),kdof2(chk))=K_c(kdof1(chk),kdof2(chk)) - 1k(chk); 
K_c(kdof2(chk),kdof1(chk))=K_c(kdof2(chk),kdof1(chk)) - 1k(chk); 
K_c(kdof2(chk),kdof2(chk))=K_c(kdof2(chk),kdof2(chk)) + Ik(chk); 
end 


changek=input(‘Are there any other changes to the stiffness matrix? l=yes 2=no ‘); 
end 


changec=input("Would you like to change your damping matrix? l=yes 2=no ‘); 
while changec==1 
chc=chc+1; % counter for determining # changes to damping matrix. 
flag_c=1; % mapping matrix flag 
Ic(chc)=input(input value of added damper (Ibs-sec/in) ‘); 
cdof1(chc)=input(‘input constrained dof where 1st end of added damper is applied ‘); 
cdof2(chc)=input(input dof where 2nd end of added damper is applied or O for 
ground '); 
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% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 


oO 
if cdof1(chc)~=cset 
cset=[cset cdof1(chc)]; 
end 
if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 


% Classical method of reformulatting the [C] matrix. 
Jo 
C_c(cdof1(chc),cdof1(chc))=C_c(cdof1(chc),cdof1(chc)) + Ic(chc); 
if cdof2(chc)~=0 
C_c(cdof1(chc),cdof2(chc))=C_c(cdof1(chc),cdof2(chc)) - Ic(chc); 
C_c(cdof2(chc),cdofi(chc))=C_c(cdof2(chc),cdof1l(chc)) - Ic(chc); 
C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc); 
end 
changec=input(‘Are there any other changes to the damping matrix? l=yes 2=no '); 
end 
changem=input(‘Would you like to change your mass matrix? l=yes 2=no '); 
while changem==1 
chm=chm+1; % counter for determining # changes to mass matrix. 
Im(chm)=input(‘input value of added mass (Ibs-sec’2/in) '); 
mdof(chm)=input(‘input constrained dof where mass is added _'); 


% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 


oO 
if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 
% Classical method of reformulatting the [M] matrix. 
% 
M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + Im(chm); 


changem=input(‘Are there any other changes to the mass matrix? l=yes 2=no ‘); 
end 


% Handling of spring added to ground 
J 


Oo 
ground = find(cset==0); 
cset(ground) = []; 
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% Formulation of iset matrix which is the set of dofs other than 

% those in the cset matrix where response information is desired. 

% 

iset=input('What dofs other than those where change occured, are you interested in? '); 


% Formation of eset matrix which is i U c. 
J 
eset=[iset cset]; 


% Choosing the FRF interested in 


Oo 
resp_dof=input(‘What response dof are you interested in? '); 
inp_dof=input(‘What input dof are you interested in? '); 
frf_index=find(eset==resp_dof | eset==inp_dof); 
if resp_dof==inp_dof 
frf_index=[frf_index frf_index]; 
end 
% Structural Synthesis method of calculating FRF 


Oo PDP ADP ABP ADP ADP DP PG AGOGO GF AGLI OLLI ALI ALL I ALLS ALLEL LILLIE FLEAS ALLE ALO. AO AGP SOG LILIA LG Ca CAO Irena) (aise 


% Calculation of unchanged FRF 

% 

freq=.01:.01:10; 

omega=2*pi* freq; 

j=length(omega); 

for w=1:} 
Z=K-+i*omega(w)*C-omega(w)*2*M; 
H=inv(Z); 


% Rearrangement of original [H] to form [Hee], [Hec], [Hec], [Hce] 
% Also provides for if user makes no changes to the model 

% 

Hee=H(eset,eset); 


if cset==[] 
Hec=0; Hcc=0; Hce=0; 
else 
Hec=Heset,cset); 
Hcc=H(cset,cset); 
Hce=H(cset,eset); 
end 


% Formation of [del_K], [del_C], [del_M], & [del_Z] 


% 
del_K=zeros(length(cset)); % intializes the change matrices 
del_C=zeros(length(cset)); % to zero and sets their sizes 
del_M=zeros(length(cset)); % as [cset x cset] matrix 
if chk~=0 

for a=1:chk 


indexk1=find(kdof1(a)==cset); 9% finds where in cset matrix 
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indexk2=find(kdof2(a)==cset); 9% the change in k occurs 


del_K(indexk1 ,indexk1)=del_K(indexk1 ,indexk1) + Ik(a); 
del_K(indexk1,indexk2)=del_K(ndexk1 ,indexk2) - lk(a); 
del_K(indexk2,indexk1)=del_K(indexk2,indexk1) - lk(a); 
del_K(indexk2,indexk2)=del_K(indexk2,indexk2) + lk(a); 
end 
end 


if chc~=0 
for a=1:chc 
indexc1=find(cdof1(a)==cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); 9% the change in k occurs 


del_C(indexc1,indexc1)=del_C(indexc1,indexc1) + Ic(a); 
del_C(indexc1 ,indexc2)=del_C(indexc1 ,indexc2) - Ic(a); 
del_C(indexc2,indexc1)=del_C(indexc2,indexc1) - Ic(a); 
del_C(indexc2,indexc2)=del_C(indexc2,indexc2) + Ic(a); 
end 
end 


if chm~=0 
for a=1:chm 
indexm=find(mdof(a)==cset); % finds where in cset matrix 
% the change in k occurs 


del_M(indexm,indexm)=del_M(indexm,indexm) + Im(a); 
end 


end 
if cset==[] 
del_Z=0; 
else 


del_Z=del_K+i*omega(w)*del_C-omega(w)2*del_M; 
end 


H_star=Hee-Hec*del_Z*inv(eye(size(Hcc,1)) + Hcc*del_Z)*Hce; 
% Capturing the FRF interested in at each frequency swept thru 
% 

hstar(w)=H_star(frf_index(1),frf_index(2)); 


% Classical method of calculating the FRF by reassembling [Z] and taking the 
% inverse. 
Jo 


Z_c=K_c+i*omega(w)*C_c-omega(w)42*M_c; 


H_c=inv(Z_c); 
H_c=H_c(eset,eset); 
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% Capturing the FRF interested in at each frequency swept thru 
% 
h_c(w)=H_c(frf_index(1),frf_index(2)); 


end 


% Plotting of FRF using Classical & Frequency Synthesis Methods 
% 
semilogy(freg,abs(h_c),'r--’,freq,abs(hstar),'b’) 
title({‘Frequency Response Function (H’,int2str(resp_dof),int2str(inp_dof),')' ]) 
ylabel(FRF Amplitude (in/ib)’) 
xlabel(‘Frequency (Hz)') 
legend(‘Classical’,’Synthesis’) 


change_model=input(‘Would you like to change your model again? l=yes 2=no '); 
end 


Ids 


% This program is designed to calculate FRFs by inverting the impedance matrix. 

% Changing the mass, stiffness, and damping matrices, new FRFs will be calculated by 
% reformulating the matrices, and using the synthesis method in the frequency domain. 
% 


clear 
% loading of the baseline structure 
oO 
disp(‘Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
[M_K_C,p]=uigetfile('*.mat','Load [M], [K], [C]’); 
load (M_K_C) 
dofs=1:ndof; 


change_model=1; 
while change_model==1; 


changek=2; changec=2; changem=2; % initializes logic variables to ‘no’ status 
flag _k=0; flag_c=0; flag _m=0; % initializes change flag 

chk=0; chc=0; chm=0; % initializes counter for # of changes to matrices 

cset=0; % initializes cset matrix to zero to avoid null index error in 


% handling the spring to ground if no changes are made to 
% the structure. 


Hey_star=[]; HV_c=[]; % (re)initializes these matrices to an empty matrix 
% to avoid matrix size differences each time the 
% model is changed. 
K_c=K; M_c=M; C_c=C; % initializes change matrices which will be 
% used in the classical method, to the original 
% matrices. 
Kb=zeros(size(K, 1)); % initializes the spring, and damper to base matrices 
Cb=zeros(size(C, 1)); % to zero. 


% Designation where base excitation is located, and the spring constant which 
: connects the substructure to the base which is moving. 
oO 
bset=input('‘At what dof(s) are your structure excited? '); 
for b_spring=1:length(bset); 
kb(b_spring)=input(fprintf(‘input the value of the spring constant which connects dof 
Jog to the base ',bset(b_spring))); 
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end 
cb=zeros(1,length(bset)); Zo initializes the damper constant which connects the 
%o substructure to the base which is moving, to zero. 


% Correcting the classical [K] & [Kb] matrices to include the spring element 
% connecting the substructure to the base(s). 


oO 
for b_spring=1:length(bset); 
K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+... 
kb(b_spring); 
Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring)) +... 
kb(b_spring); 
end 


changek=input(‘Would you like to change your stiffness matrix? 1=yes 2=no '); 
while changek==1 
chk=chk+1; % counter for determining # change of K matrix 
flag_k=1; 
1k(chk)=input(‘input value of added spring (lbs/in) '); 
kdof1(chk)=input(‘input constrained dof where 1st end of added spring is applied '); 
kdof2(chk)=input(‘input dof where 2nd end of added spring is applied, b for base, or 
0 for ground ‘); 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 
Jo 
if kdof1(chk)~=bset & kdof2(chk)=='b' 
while kdof1(chk)~=bset & kdof2(chk)=="b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(If this is what you desire, go back and include this dof as an excitation 
point.') 
model_change=input(Was this change correct? I1=yes 2=no ’); 


if model_change== 
error('Go back and change your basic model.’) 

else 
kdof1(chk)=input(‘input constrained dof where 1st end of added damper is 

applied ’); 
kdof2(chk)=input(‘input dof where 2nd end of added damper is applied, b 
for base, or O for ground ‘); 
end 
end 
end 


% updates the base spring constant and does not count these changes 
% in the cset matrix. 


TiS 


% 

if (find(kdof1(chk)==bset))~=[] & kdof2(chk)=='b' 
base=find(kdof1(chk)==bset); 
kb(base) = kb(base) + Ik(chk); 


se 
if kdof1(chk)~=cset % comparison of change dof with the c-set 
cset=[cset kdof1(chk)]; % matrix and formation of new c-set matrix. 
end 
if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 
end 


% Classical method of reformulatting the [K] matrix. 


O 
K_c(kdof1(chk),kdof1(chk))=K_c(kdof1(chk),kdof1(chk)) + lk(chk); 
if kdof2(chk)~=0 & kdof2(chk)~="b' 
K_c(kdof1(chk),kdof2(chk))=K_c(kdof1(chk),kdof2(chk)) - 1k(chk); 
K_c(kdof2(chk),kdof1(chk))=K_c(kdof2(chk),kdof1 (chk)) - 1k(chk); 
K_c(kdof2(chk),kdof2(chk))=K_c(kdof2(chk),kdof2(chk)) + lk(chk); 
end 


% Classical method of reformulatting [Kb] matrix. 
%o 

if kdof2(chk)=="b' 

kbb=lk(chk); 

Kb(kdofl (chk), kdof1(chk))=Kb(kdof1 (chk),kdof1(chk)) + kbb; 
end 


changek=input(‘Are there any other changes to the stiffness matrix? 1=yes 2=no '); 
end 


changec=input(‘Would you like to change your damping matrix? 1=yes 2=no '); 
while changec==1 
chc=che+1; % counter for determining # changes to damping matrix. 
flag_c=1; 
lc(chc)=input(‘input value of added damper (Ibs-sec/in) '); 


cdof1(chc)=input(‘input constrained dof where 1st end of added damper is applied _‘); 
cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b for base, or 
0 for ground ‘); 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 


O 
if cdof1(chc)~=bset & cdof2(chc)=="b' 
while cdof1l(chc)~=bset & cdof2(chc)=="b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(‘If this is what you desire, go back and include this dof as an excitation 
point.’) 
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model_change=input('Was this change correct? 1=yes 2=no '); 
if model_change==1 


error(‘Go back and change your basic model.) 
else 


cdof1(chc)=input(‘input constrained dof where 1st end of added damper is 


applied °); 
cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground '); 
end 
end 
end 


% Updates the base damper constant and does not count these changes 
% in the cset matrix. 


% 
if (find(cdof1(chc)==bset))~=[] & cdof2(chc)=="b' 
base=find(cdof1(chc)==bset); 


cb(base) = cb(base) + Ic(chc); 
else 


if cdof1(chc)~=cset 


% comparison of change dof with the c-set 
cset=[cset cdof1(chc)]; 


% matrix and formation of new c-set matrix. 
end 
if cdof2(chc)~=cset 


cset=[cset cdof2(chc)]; 
end 


end 


% Classical method of reformulatting the [C] matrix. 

% 

C_c(cdof1(chc),cdof1(chc))=C_c(cdofl(chc),cdof1(chc)) + lc(chc); 
if cdof2(chc)~=0 & cdof2(chc)~="b' 


C_c(cdof1(chc),cdof2(chc))=C_c(cdof1(chc),cdof2(chc)) - le(chc); 

C_c(cdof2(chc),cdof1(chc))=C_c(cdof2(chc),cdof1(chc)) - lc(chc); 

C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc) 
end 


> 


% Classical method of reformulatting [Cb] matrix. 
% 


oO 
if cdof2(chc)=="b' 
cbb=lc(chc); 


Cb(cdof1(chc),cdof1(chc))=Cb(cdof1(chc),cdof1(chc)) + cbb; 
end 
changec=input(‘Are there any other changes to the damping matrix? l=yes 2=no ’); 
end 


changem=input('Would you like to change your mass matrix? 1=yes 2=no ‘); 
while changem==1 


chm=chm+1; 
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flag_m=1; 
lm(chm)=input(‘input value of added mass (lbs-sec’2/in) '); 
mdof(chm)=input(‘input constrained dof where mass is added '); 


% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 


O 
if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 


% Classical method of reformulatting the [M] matrix. 
% 
M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + lm(chm); 


changem=input(‘Are there any other changes to the mass matrix? l=yes 2=no '); 
end 


% Handling of spring added to ground 
J 


‘0 
ground = find(cset==0); 
cset(ground) = []; 


% Formulation of iset matrix which is the set of dofs other than 

% those in the cset matrix where response information is desired. 

% 

iset=input("What dofs other than those where change occured, are you interested in? ‘); 


% Formation of zset (c U b) and eset (1 U c U b) matrices. 
% 
zset=[cset bset]; eset=[iset cset bset]; 


% Choosing the FRF interested in 

%o 

resp_dof=input(‘What response dof are you interested in? ‘); 
frf_index=find(eset==resp_dof); 


% Choosing the type of input and out FRF interested in 
O 
inp_type=input(Is the base excitation in the form of: 1=displacement or 2=acceleration 
resp_type=input("What type of response are you interested in: 1=displacement or 
2=acceleration '); 


% Structural Synthesis method of calculating FRF 


% 
% Calculation of unchanged FRF 
% 
freq=.01:1:250+.01; 
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omega=2*pi*freq; 

j=length(omega); 

for w=1:] 
Z=K+i*omega(w)*C-omega(w)2*M; 
H=inv(Z); 


% Rearrangement of original [H] to form [Heb], [Hez], [Hzz], [Hzb] 
Ye 


oO 
Heb=H(eset,bset); 
Hez=H(eset,zset); 
Hzz=H(zset,zset); 
Hzb=H(zset,bset); 


% Formation of [del_Kc], [del_Cc], [del_Mc], [del_Zc], [del_Zb] 
% & [del_Z} 
Te 


oO 
del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Cc=zeros(length(cset)); % to zero and sets their sizes 
del_Mc=zeros(length(cset)); % as [cset x cset] matrix 
if flag_k==1 
for a=1:chk 
if (find(kdof1(a)==bset))~=[] & kdof2(a)=="b' 
indexk 1=[]; % provides for base changes to not be 
indexk2=[]; % included in del_Kc matrix 
else 
indexk 1=find(kdof1(a)==cset); % finds where in cset matrix 
indexk2=find(kdof2(a)==cset); % the change in k occurs 
end 


del_Kc(ndexk1,indexk1)=del_Kc(indexk1,indexk1) + lk(a); 
del_Kc(indexk1,indexk2)=del_Kc(indexk1,indexk2) - Ik(a); 
del_Kc(indexk2,indexk1)=del_Kc(indexk2,indexk1) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 


end 
end 
if flag_c==1 
for a=1:che 
if (find(cdof1(a)==bset))~=[] & cdof2(a)=="b' 
indexc1=[]; 
indexc2=[]; 
else 
indexc1=find(cdof1(a)==cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); % the change in C occurs 
end 


del_Cc(indexc1,indexcl1 )=del_Cc(indexc1,indexc1) + Ic(a); 
del_Cc(indexc1,indexc2)=del_Cc(indexc1,indexc2) - 1c(a); 
del_Cc(indexc2,indexc1)=del_Cc(indexc2,1indexc1) - Ic(a); 
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del_Cc(indexc2,indexc2)=del_Cc(indexc2,indexc2) + lc(a); 


end 

end 

if flag_m==1 
for a=1:chm 


indexm=find(mdof(a)==cset); 9% finds where in cset matrix 
% the change in k occurs 


del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + lm(a); 
end 
end 


del_Zb=diag(kb+1*omega(w)*cb); % diag is used in the event >1 base 
% excitation is present 


del_Zc=del_Kc-*omega(w)*del_Cc-omega(w)2*del_Mc; 


del_Z=[ del_Zc zeros(size(del_Zc,1),size(del_Zb,2)) 
zeros(size(del_Zb,1),size(del_Zc,2)) del_Zb]; 


% Base excitation Frequency synthesis equation 


Oo 
Hey_star(:,w)=(Heb*del_Zb -Hez*del_Z*inv(eye(size(Hzz, 1))+Hzz*del_Z)*... 
Hzb*del_Zb)*ones(1,length(bset))’; 


if inp_type==1 & resp_type==2 

Hey_star(:,w) = Hey_star(:,w)*(-omega(w)2); 
elseif inp_type==2 & resp_type==1 

Hey_star(:,w) = Hey_star(:,w)*(-1/omega(w)2); 
end 


% Classical method of “gpeuae the FRF by ee [Z] and taking the 
% inverse. 
To 


Il=eye(size(M_c,1)); 

=ones(1,ndof)’; 
Z_c=(K_c+*omega(w)*C_c-omega(w)42*M_c); 
H_c=inv(Z_c); 
HV_c(:,w) = (H_c*(Kb+i*omega(w)*Cb))* YV; 


if inp_type==1 & resp_type==2 
HV_c(:,w) = HV_c(:,w)*(-omega(w)2); 
elseif 1 ee type==2 & resp_type==1 
HV_c(:,w) = HV_c(:,w)*(-1/omega(w)2); 
end 
end 
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% Plotting of FRF using Classical & Frequency Synthesis Methods 


Jo 


semilogy(freq,abs(HV_c(resp_dof,:)),'r--',freq,abs(Hey_star(frf_index,:)),'b’) 
title([Frequency Response Function (H’,int2str(resp_dof),')' }) 
if inp_type==1 & resp_type==2 
ylabel(FRF Amplitude (g/in)') 
elseif inp_type==2 & resp_type==1 
ylabel(FRF Amplitude (in/g)’) 


else 
ylabel(FRF Amplitude (unitless)’') 
end 
xlabel(‘Frequency (Hz)') 


legend(‘Classical’,'Synthesis') 


change_model=input(‘Would you like to change your model again? l=yes 2=no ‘); 
end 
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fsynmain.m 


% This program is designed to calculate the new FREs of a system after changes in the 
% mass, stiffness, and damping matrices have been made. The new FRFs will be 
% calculated using the synthesis method in the frequency domain. 


%o 
clear 


% loads the baseline structure to be modified, allows for modifications to the baseline 
% structure, and calls the function delta.m 


Jo 

load change 

% Designation of the frequency range and step size 

qe initial stepsize end 

Coes =[{ 01 on — 250+.01 ]; 


freq=freq_input(1):freq_input(2):freq_input(3); 
omega=2*pi*freq; 


[wn,phi,zeta,Mmodal,phi_norm,Cmodal] = modal(M,K,C); 
[hee] = frfmodal(wn,phi,zeta,Mmodal,omega,eset); 


[h_star] = frfsynth(hee,kb,cb,omega,excite,eset,iset, bset,del_Kc,del_Cc,del_Mc); 


% Choosing the type of input and output FRF interested in 
%o 
if excite==2 
inp_type=input(Ts the base excitation in the form of: 1=displacement or 2=acceleration 


resp_type=input(‘What type of response are you interested in: 1=displacement or 
2=acceleration '); 
if inp_type==1 & resp_type==2 
h_star = h_star * (-omega(w)2); 
elseif inp_type==2 & resp_type==1 
h_star = h_star * (-1/omega(w)2); 
end 
end 


resp_dof=input(‘What response dof are you interested in? ‘); 


if excite==2 

inp_dof={]; 
else 

inp_dof=input(‘What input dof are you interested in? ‘); 
end 


122 


[H_star_desired] = frf_sift(eset,resp_dof,inp_dof,excite,h_star); 


semilogy(freq',abs(H_star_desired),'g') 
title({‘Synthesized Frequency Response Function H’,int2str([resp_dof inp_dof]) J) 
if inp_type==1 & resp_type== 
ylabel(FRF Amplitude (g/in)’) 
elseif inp_type==2 & resp_type==1 
ylabel(FRF Amplitude (in/g)’) 
else 
ylabel(FRF Amplitude (unitless)’) 
end 
xlabel(‘Frequency (rad/s)') 
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% This program provides for the loading of the baseline structure, accepts the user changes 
% to the [M], [K], and [C] matrices, calls the function delta.m, and stores this information 
%o to be loaded from another program. 


% 
clear 


changek=2; changec=2; changem=2; % initializes logic variables 


% to ‘no’ 
status 
flag k=0; flag_c=0; flag_m=0; 
chk=0; chc=0; chm=0; % initializes counter for # of changes to matrices 
cset=0; bset=[]; % initializes cset matrix to zero to avoid 


% null index error if no changes are made to 
% the structure and bset to empty matrix to 
% prevent error in the event this is not a 

% base excitation and bset does not exist. 


kdofl={]; kdof2={]; % initalizes the matrices which keep 
cdofl={]; cdof2={J; % track of the dofs where changes occur 
mdof={]; 


Ik=0; Ic=0; Im=0; % initializes the value of added prarmeters 


kb0O=0; cb0=0; % initializes kb and cb to zero to prevent error in the event 
% this 1s not a base excitation and kb & cb does not exist. 


b=num2str('b'); % allows b to be entered as a numerical value, but 
% declares b a string variable to used for 
% comparisons in logic statements. 


disp(‘Are there already reduced FRF & IRF matrices in the correct format available,’) 
disp(‘or does the presynthesized FRF &IRF matrices need to be generated using the dof ') 
disp(‘locations where you desire to change the structure? ') 

FRF_TRF=input('l=reduced FRF/IRF already exists 2=generate reduced FRF/IRF ‘'); 


% Protects against user making an error in choosing how FRF & IRF is obtained. 


oO 
if FRF_IRF~=[1 2] 
while FRF_IRF~=[1 2] 
disp(Error in choosing how FREF/IRF is obtained. Choose 1 or 2.’) 
FRF_IRF=input('1=reduced FRF/IRF already exists 2=generate reduced FRF/IRF ‘); 
end 
end 
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% Loading the presynthesized FRF/IRFs, corresponing frequency vector, and vector of 
% dofs in the eset if reduced FRF/IRF already exists. Loading of M, K, & C matrices if 
% FRF/IRF needs to be generated 


oO 

if FRF_IRF==1 
disp(‘Select the file which contains your presynthesized FRF/IRFs, a corresponding 

frequency ') 
disp(‘vector, and the vector of all the dofs that will be in your eset vector.’) 
pause(2) 
[hee_dofs_omega,p]=uigetfile('*.mat',"Load [FRF], [IRF], {dofs}, {omega}'); 
load (hee_dofs_omega) 


else 
disp(‘Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
[M_K_C,p]=uigetfile('*.mat',Load [M}, [K}, [C}’); 
load (M_K_C) 
dofs=1:ndof; 
end 


excite=input(‘How is the structure excited? l=system 2=base '); 
% Protects against user making an error in choosing the type of excitation. 


oO 
if excite~=[1 2] 
while excite~=[1 2] 
disp(Error in choosing type of excitation. Choose 1 or 2.) 
excite=input(‘How is the structure excited? 1=system 2=base '); 
end 
end 


if excite==2 
% Designation where base excitation is located. 
% 
bset=input(‘At what dof(s) are your structure excited? ‘); 


% Protects against user making an error in choosing the location of the base excitation. 
Jo 
for b_dof=1:length(bset) 
if bset(b_dof)~=dofs 
while bset(b_dof)~=dofs 
fprintf(‘Error in choosing location of base excitation at dof %g.\n’.... 
bset(b_dof)) 
disp(‘Either the dof chosen does not exist on your structure, or is not ’) 
disp(‘included in your list of dofs for the eset vector. ‘) 
bset(b_dof)=input(fprintf(‘Choose again the dof for base exitation %g ".... 
b_dof)); 
end 
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end 
end 


% Designation of the spring constant which connects the substructure to the base and 
% protection against user making an error in choosing the value of the constant(s). 


oO 
for b_spring=1:length(bset) 
kbO(b_spring)=input(fprintf(‘input the value of the spring constant which connects dof 
Jog to the base ',bset(b_spring))); 
kbO_zero=find(kb0==0); 


if length(kb0)~=b_spring | kbO_zero~=[] 
while length(kb0)~=b_spring | kbO_zero~=[] 
fprintf(‘You made an error in entering the value for dof %g base 
excitation spring constant.\n',bset(b_spring)) 

kb0(b_spring)=0; 

kbO(b_spring)=input(fprintf(‘Re-input the value of the spring constant 

which connects dof %g to the base ',bset(b_spring))); 

kb0_zero=find(kb0==0); 


end 
end 
end 
cb0=zeros(1,length(bset)); % initializes the damper constant which connects the 
% substructure to the base which is moving, to zero. 
end 


changek=input(‘Would you like to make stiffness modifications to the structure? l=yes 
2=no '); 
while changek==1 


% Designation of the spring change and protection against user making an error in 
% choosing the value of the constant(s). 


0 
chk=chk+1; % counter for determining # changes to stiffness matrix. 
lk(chk)=input(‘input value of added spring (lbs/in) ‘); 
if length(lk)~=chk 
while length(k)~=chk 
disp("You made an error entering the value for the stiffness change.") 
lk(chk)=0; 
1k(chk)=input(‘Re-input value of added spring (lbs/in) '); 
end 
end 


% User selection of where stiffness changes occur. Protects against user making an 
% error in choosing the locations of spring changes. 


O 
kdof1(chk)=input(‘input dof where ist end of added spring is applied ‘); 
if kdof1(chk)~=dofs 
while kdof1(chk)~=dofs 
disp(‘Error in choosing location of 1st end of added spring. Either the dof 
chosen does’) 
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disp(‘not exist on your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ' ) 
kdof1(chk)=input(‘Re-input constrained dof where Ist end of added spring is 
applied ’); 
end 
end 


kdof2(chk)=input(‘input dof where 2nd end of added spring is applied, b for base, or 0 
for ground '); 
if kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~='b' 
while kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~="b' 
disp(Error in choosing location of 2nd end of added spring. Either the dof 
chosen does’) 
disp(not exist on your structure, or 1s not included in your list of dofs for the’) 
disp(‘eset vector. ' ) 
kdof2(chk)=input(‘Re-input dof where 2nd end of added spring is applied, b for 
base, or 0 for ground '); 
end 
end 


if excite==1 & kdof2(chk)=='b' 
while excite==1 & kdof2(chk)=="b' 
disp(Error in choosing location of 2nd end of added spring. You did not 
, designate this’) 
disp(‘as a base excitation structure. Therefore b is not a dof option.") 
kdof2(chk)=input('Re-input dof where 2nd end of added spring is applied, b for 
base, or 0 for ground ‘); 
end 
end 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should be handled when asked where the structure is excited. 
% 
if excite==2 & kdof1(chk)~=bset & kdof2(chk)=="b' 
while kdof1(chk)~=bset & kdof2(chk)=="b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(If this is what you desire, go back and include this dof as an excitation 
point.') 
model_change=input(‘Was this change correct? 1=yes 2=no '); 
if model_change==1 
error(‘Go back and change your basic model.’) 
else 
kdof] (chk)=input(‘Re-input constrained dof where 1st end of added spring is 
applied ’); 
kdof2(chk)=input(‘Re-input dof where 2nd end of added spring is applied, b 
for base, or 0 for ground '); 
end 
end 
end 
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% Provision for base spring constant changes to not be included in the cset matrix. 


% 
if excite==2 & (find(kdof1(chk)==bset))~=[] & kdof2(chk)=="'b' 
cset=cset; 
else 
if kdof1(chk)~=cset % comparison of change dof with the c-set 
cset=[cset kdof1(chk)]; % matrix and formation of new c-set matrix. 
end 
if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 
end 


changek=input(‘Are there any other stiffness modifications to the structure? 1=yes 2=no’); 
end 


changec=input(‘Would you like to make damping modifications to the structure? 1=yes 
2=no '); 
while changec==1 


% Designation of the damper change and protection against user making an error in 
% choosing the value of the constant(s). 
% 
che=chce+]; % counter for determining # changes to damping matrix. 
lc(chc)=input(‘input value of added damper (Ibs-sec/in) ‘); 
if length(Ic)~=che 
while length(c)~=chc 
disp('You made an error entering the value for the damper change.’) 
Ic(chc)=0; 
lc(chc)=input('Re-input value of added damper (lbs-sec/in) '); 
end 
end 


% User selection of where damper changes occur. Protects against user making an error 
% in choosing the locations of damper changes. 


b 
cdof1(chc)=input(‘input constrained dof where 1st end of added damper is applied ‘); 
if cdof1 (chc)~=dofs 
while cdof1(chc)~=dofs 
disp(Error in choosing location of 1st end of added damper. Either the dof 
chosen does’) 
disp(‘not exist on your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ‘ ) 
cdof1(chc)=input(‘Re-input constrained dof where 1st end of added damper is 
applied ’); 
end 
end 


cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b for base, or 0 
for ground ‘); 
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if cdof2(chc)~=dofs & cdof2(chc)~=0 & cdof2(chc)~="b' 
while cdof2(chc)~=dofs & cdof2(chc)~=0 & cdof2(chc)~='b' 
disp(Error in choosing location of 2nd end of added damper. Either the dof 
chosen does’) 
disp(not exist On your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ' 
cdof2(chc)=input(‘Re-input dof where 2nd end of added damper is applied, b for 
base, or QO for ground '); 
end 
end 


if excite==1 & cdof2(chc)=="b' 
while excite==1 & cdof2(chc)=="b' 
disp(Error in choosing location of 2nd end of added damper. You did not 
designate this’) 
disp(as a base excitation structure. Therefore b is not a dof option.') 
cdof2(chc)=input('Re-input dof where 2nd end of added damper is applied, b for 
base, or O for ground '); 
end 
end 


% Let's the user know that he 1s not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from these dof to base constitutes an addition 
% of excitation points. This should be handled when asked where the structure is 
% excited. 
%o 
if excite==2 & cdofl(chc)~=bset & cdof2(chc)=='b' 
while cdof1(chc)~=bset & cdof2(chc)=="b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(‘If this is what you desire, go back and include this dof as an excitation 
point.') 
model_change=input('Was this change correct? 1=yes 2=no '); 


if model_change== 
error('Go back and change your basic model.’) 

else : 
cdof1(chc)=input(‘Re-input constrained dof where 1st end of added damper is 

applied ’); 
cdof2(chc)=input(‘Re-input dof where 2nd end of added damper is applied, b 
for base, or OQ for ground ');_ 
end 
end 
end 


% Provision for base damper constant changes to not be included in the cset matrix. 


O 
if excite==2 & (find(cdof1(chc)==bset))~=[] & cdof2(chc)==b' 


cset=cset; 
else 
if cdof1(chc)~=cset % comparison of change dof with the c-set 
cset=[cset cdof1(chc)]; % matrix and formation of new c-set matrix. 


129 


end 
if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 


changec=input(‘Are there any other damping modifications to the structure? l=yes 2=no’); 
end 


changem=input("‘Would you like to make a mass modification to the structure? l=yes 
2=no ); 
while changem—1 


% Designation of the mass change and protection against user making an error in 
choosing the 

% value of the constant(s). 

Jo 

chm=chm-+1; % counter for determining # changes to mass matrix. 

Im(chm)=input(‘input value of added mass (Ibs-sec’2/in) '); 

if length(lm)~=chm 

while length(lm)~=chm 

disp(‘You made an error entering the value for the mass change.) 
Im(chm)=0; 
Im(chm)=input(‘Re-input value of added mass (lbs-sec’2/in) '); 
end 
end 


% User selection of where mass changes occur. Protects against user making an error in 
% choosing the locations of mass changes. 


J 
mdof(chm)=input(‘input constrained dof where mass is added '); 
if mdof(chm)~=dofs 
while mdof(chm)~=dofs 
disp(Error in choosing location of mass change. Either the dof chosen does 
not’ 
disp(‘exist on your structure, or is not included in your list of dofs for the’) 
disp(eset vector. ' ) 
mdof(chm)=input(‘input constrained dof where mass is added ‘); 
end 
end 
if mdof(chm)~=cset % comparison of change dof with the c-set matrix 
cset=[cset mdof(chm)]; % and formation of new c-set matrix. 
end 


changem=input('Are there any other mass modifications to the structure? l=yes 2=no ’); 
end 


% Handling of spring added to ground to eliminate zero from cset 
oO 
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ground = find(cset==0); 
cset(ground) = []; 


disp(‘The following is a list of dofs were you have made changes to your structure.’) 
disp(‘This represents your {cset} vector:') 
disp(cset) 


% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 
%o 
if FRF_IRF==1 

disp(‘Since you inputed your own FRF/IRF matrices, your iset is chosen as all dofs 
remaining ') 

disp(‘in your eset/dof vector after extracting the cset vector.’) 

iset=dofs; 

for a=1 :length(cset) 

extract(a)=find(cset(a)==dofs); 

end 

iset(extract)=[]; 
else 

iset=input("What dofs other than those where change occured, are you interested in? '); 
end 


% Formation of zset (c U b) and eset (1 U c U b) matrices. 
% 


zset=[cset bset]; eset=[iset cset bset]; 


{del_Kc,del_Cc,del_Mc,kb,cb] = 
delta(cset, bset,kdof1 ,kdof2,]k,cdof1 ,cdof2,lc,mdof,lm,kb0,cb0); 


save change M K C ndof del_Kc del_Cc del_Mc kb cb iset cset bset zset eset excite 
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function [wn,phi,zeta,Mmodal,phi_norm,Cmodal] = modal(M,K,C) 


% This function is designed to calculate the natural frequencies and 
% mode shapes 
% 


oO 

ndof=size(M, 1); 
[PHI,lam]=eig(M\K); 
for j=1:ndof; 

lambda(j)=lam(j,}); 
end 
[wn,I]=sort(sqrt(lambda)); 
wn=rot90(wn,-1); 
for j=1:ndof; 

phi(:,j)=PHIC:,1Q)); 
en 


% Formulation of modal mass and damping matrices 


O 
Mmodal=phi'*M*phi; 
phi_norm=phi/(Mmodal“.5); 
Mmodal=diag(Mmodal); 
Cmodal=phi'*C* phi; 
Cmodal=diag(Cmodal); 
%ozeta0=0.0; 
%zeta=zeta0*ones(1,length(Cmodal)); % Constant modal damping ratio 


zeta=Cmodal./(2*Mmodal.*wn); Jo Modal damping ratio derived 
% from proportional damping matrix [C] 
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frfmodal.m 


function [h] = frfmodal(wn,phi,zeta,Mmodal,omega,frf_dof) 
% This function is designed to calculate the FRF using mode shapes and 


% summing the FRF over a number of modes. 
% 


ncol=(length(frf_dof)*(length(frf_dof)+1))/2; 


% Calculation of frequency response function 


% 
h=zeros(length(omega),ncol); % Initialization to zero, and sizing 
% of the freq x lower 

triangle matrix. 
h_old=0; h_new=0; % initialization of the matrices to 

% be used to determine when modal 

% convergence has occured. 
h_check=1; % Initialization of the stop criteria for summing the 


% modes 


nmodes=size(phi); b=0; % Defining the number of modes that exist and 
% initializing b which keeps track of the 
% of modes used for convergence of the FRF. 


HOdprime_aa=zeros(length(frf_dof)); 
for b=1:nmodes 


% Elimination of unnessecary elements and rearrangement of original 

% mode shapes {phi} to form {phi_red} and [num_red] which is now an 

% frf_dof x frf_dof matrix. 

J 

num_red=phi(frf_dof,b)*phi(frf_dof,b)’; 

for w=1:length(omega) 
den=(wn(b)42 - omega(w)2 + 2*i*zeta(b)* wn(b)*omega(w))*Mmodal(b); 
H_red=num_red./den; 


% Changes the [H_red] FRFs from an frf_dof x frf_dof matrix to 
% a freq x lower tnangle matrix. 


Jo 
H_lowtmi=til(H_red); % pulls out the lower triangle 

% and places zeros everywhere else 
symtry = find(H_lowtm==0); % finds zero values in the 
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H_lowtri(symtry) = []; 


h_mode(w,:)=H_lowt, 


end 
h =h + h_mode; 


end 


% lower triangle matrix 


% deletes all values of zero and 

% turns the remaining elements 

% into a 1 x length(lower triangle) 
% vector 


% saves the lower triangle vector 
% at each freq. These FRFs are for 
% a single mode only. 


% Summing of each modal FRF 
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frfsynth.m 


function [h_star] = frfsynth(hee,kb,cb,omega,excite,eset,iset, bset.... 
del_Kc,del_Cc,del_Mc) 


% This function is designed to calculate the new FRF using the synthesis 
% method. 
% 


for w=1:length(omega) 
count=0; % initializes the counter which keeps track of 
% which column of the hee matrix is to be 
% placed into the refomation of matrix [Hee] 


for col = 1:length(eset) % reforms the [Hee] lower 
for row = col:length(eset) % triangle matrix 
count=count + 1; 
Hee_lowtri(row,col)=hee(w,count); 
end 
end 


% reforms the original symmetric [Hee] matrix 

%o 

Hee=Hee_lowtn; 

[sym_row,sym_col] = find(Hee_lowtri==0); 

for n=1:length(sym_row) 
Hee(sym_row(n),sym_col(n))=Hee(sym_col(n),sym_row(n)); 

end 


e=1:length(eset); 
c=length(iset)+1:length(eset)-length(bset); 
b=length(eset)-length(bset)+1 :length(eset); 
z=length(iset)+1 :length(eset); 


% Rearrangement of [Hee] to form [Hec], [Hcc], [Hce] if the structure 
% is system excited. 
Jo 
if excite==1 
Hec=Hee(e,c); 
Hce=Hee(c,e); 
Hcc=Hee(c,c); 
end 


% Rearrangement of [Hee] to form [Heb], [Hez], [Hzz], [Hzb] if the 
% structure is base excited. 
%o 
if excite== 
Heb=Hee(e,b); 
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Hez=Hee(e,z); 

Hzz=Hee(z,z); 

Hzb=Hee(z,b); 
end 


% Formation of [del_Zc], [del_Zb] & [del_Z] 
%o 
del_Zc=del_Kc+i*omega(w)*del_Cc-omega(w)42*del_Mc; 


if excite==2 
del_Zb=diag(kb+i*omega(w)*cb); % diag is used in the event 
% >1 base excitation 
% is present 
del_Z=[ del_Zc zeros(size(del_Zc,1),size(del_Zb,2)) 


zeros(size(del_Zb, 1),size(del_Zc,2)) del_Zb]; 
end 


% Determining which synthesis equations are to be used 
Jo 
if excite== 


% Structure excitation Frequency synthesis equation 
0 
H_star=Hee-Hec*del_Zc*inv(eye(size(Hcc,1)) + Hcc*del_Zc)*Hce; 
else 


% Base excitation Frequency synthesis equation 
% 


h_star(:,w)=(Heb*del_Zb - 
Hez*del_Z*inv(eye(size(Hzz, 1))+Hzz*del_Z)*Hzb*del_Zb)*ones(1 ,length(bset))’; 
end 
if excite==1 
% Changes the [H*] FRFs from an eset x eset matrix to a freq x lower 
% triangle matrix, and creates an index matrix which keeps track of 
% the original coordinates of the FRFs in the [H*] matrix. 


O 


Hstar_lowtn=tril(H_star); % pulls out the lower triangle and places 
% zeros everywhere else 


symtry = find(Hstar_lowtri==0); % finds zero values in the lower 
% triangle matrix 


Hstar_lowtri(symtry) = []; % deletes all values of zero and turns the 
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% remaining elements into a 1 x length(lower 
% triangle) vector 


h_star(w,:)=Hstar_lowtri; % saves the lower triangle vector at each freq 


end 
end 


if excite==1 
% changes h_star matrix so freqs are in each column and the FRFs are 
% in the rows 
% 
h_star=flipud(rot90(h_star)); 


end 
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function [H_star_desired] = frf_sift(eset,resp_dof,inp_dof,excite,h_star) 


% Locates the positions in the eset that the response and input dofs are refernng. 
% These locations correspond to the matrix location row and column of the desired 
% FRF in the [H*] matrix which is eset x eset and partitioned. Or these locations 
% refer to the row of the [H*] matrix which is eset x 1 for a base excitation. 

% These indices sift through the [h*] matrix rows to find the desired FRF. 

Jo 


% Creates an index matrix which keeps track of the original coordinates of the FRFs 
% in the [H_star] matrix. 


Jo 
count=0; % initializes the counter which keeps track of 
% which column of the h matrix the lower tnangle [H_red] 
J FRF went into 
for col = 1:length(eset) % creates a matrix which matches up 
for row = col:length(eset) % each element of the lower triangle 
count=count + 1; % H_star matrix with the rows in the 
hstar_index(count,:)=[row col]; % h_star matrix that they were placed in. 
end 
end 


frf_index=find(eset==resp_dof | eset==inp_dof); 
Yo prevents index error in the event response dof and input dof are the same 
%o 
if resp_dof==inp_dof 
frf_index=[frf_index frf_index]; 
end 
if excite==1 


% Uses hstar_index and finds the row in the lower triangle x freq matrix that the 
% desired FRF is located. 
qe 


oO 
frf_row=find((frf_index(1)==hstar_index(:,1) & frf_index(2)==hstar_index(:,2)) | 
(frf_index(1)==hstar_index(:,2) & frf_index(2)==hstar_index(:,1))); 


H_star_desired = h_star(frf_row,:); 


else 
H_star_desired = h_star(frf_index,:); 


end 
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APPENDIX B. STATIC DISPLACEMENT SYNTHESIS COMPUTER 


CODES 


The following is a list and a brief description of the main MATLAB 


computer codes that were written in order to perform the static displacement 


synthesis calculations. 


Guyan _ xstat.m - uses Guyan reduction methods to calculate the static 
displacements on a structure. 

ssynmain.m - the main program which calls the modularized static 
displacement synthesis programs and does the post processing of the 
results. 

statchange.m - loads the baseline structure to be modified, allows for 
modifications to the baseline structure, and calls the function statdelta.m 
statdelta.m - forms the modification matrices which will be used to form 
impedance matrices. 

modal.m - solves for the structure’s natural frequencies, mode shapes and 
other modal matrices and vectors. 

frfmodal.m - calculates the structure’s FRFs. 

statsynth.m - performs the synthesis on the baseline structure and returns 


the synthesized static displacements. 


The full codes are contained on subsequent pages. Any codes that were 
previously presented will not be repeated. 
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Guyan _ xstat.m 


% The purpose in this program is to use Guyan reduction methods in order to calculate the 
% static sisplacements on a structure. 


% 

clear 

load fe_add 

% Stiffness changes to the structure 

Jo 

K(3,3)=K(3,3)- 1000; K(18,18)=K(18,18)- 1000; 
K(93,93)=K(93,93)-1000; K(108,108)=K(108,108)-1000; 
K(3,3)=K(3,3)+1e4; K(18,18)=K(18,18)+1¢4; 
K(93,93)=K(93,93)+1e4; K(108,108)=K(108,108)+1e4; 
K(66,66)=K(66,66)+500; K(109,109)=K(109,109)+500; 
en 109)=K(66, 109)-500; K(109,66)=K(109,66)-500; 


--oneeeiee 1),1)*(-386.4); 

for rot=1:3:size(M,1)-3 

g(rot)=0; 

end 

for rot=2:3:size(M,1)-2 
g(rot)=0; 

end 


F=M*g; 
dofset=1:ndof; 


oset=dofset; 

aset=[3 18 93 108 66 109]; 

for a=1:length(aset) 
kill(a)=find(aset(a)==dofset); 

end 

oset(kill)=[]; 


nset=[aset oset]; 


Knn=K(nset,nset); | Fnn=F(nset); 

Kaa=K (aset,aset); Kao=K (aset,oset); Koo=K(oset,oset); Koa=K(oset,aset); 
Toa=(-1)*inv(Koo)* Koa; 

T=[eye(length(aset));Toal]; 


Kred_Guyan=T'*Knn?T; Fred_Guyan=T'*Fnn; 
xstat_Guyan=Kred_Guyan\Fred_Guyan; 
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ssynmain.m 


% This program is designed to calculate the new static displacements of a system after 
% changes in the mass, stiffness, and damping matrices have been made. The new static 
% displacements will be calculated using the static displacement synthesis method 

% 

clear 


load statchange 
stat_omega=. 1:.1:.4; 
C=zeros(size(C)); 


[stat_wn,stat_phi,stat_zeta,stat_Mmodal] = modal(M,K,QC); 
tO = clock; 
[stat_hee] = frfmodal(stat_wn,stat_phi,stat_zeta,stat_Mmodal,stat_omega,eset); 


[z_stat,Fred,Kred,Mred,Hstar0_dp] = statsynth(stat_hee,kb,stat_omega,eset,iset, bset,... 
del_Kc,del_Mc); 

tl = clock; 

stat_syn_time=etime(t1 ,t0) 
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statchange.m 


% This program provides for the loading of the baseline structure, accepts the user changes 
% to the [M], [K], and [C] matrices, calls the function statdelta.m, and stores this 

% information to be loaded from another program. 

% 


clear 


changek=2; changec=2; changem=2; ~% initializes logic variables ' 
% to 'no' 


status 
flag k=0; flag_c=0; flag_m=0; 
chk=0; chc=0; chm=0; % initializes counter for # of changes 
% to matrices 

cset=0; bset=[]; % initializes cset matrix to zero to avoid 

% null index error if no changes are made to 

% the structure and bset to empty matrix to 

% prevent error in the event this is not a 

% base excitation and bset does not exist. 
kdofl=[]; kdof2=[]; % initalizes the matrices which keep 
cdofl=f}; cdof2=f}; % track of the dofs where changes occur 
mdof=[]; 
lk=0; Ic=0; Im=0; % initializes the value of added prarmeters 
kbO=0; cb=0; % initializes kb and cb to zero to 

% prevent error in the event this is not a 

% base excitation and kb & cb does not exist. 
b=num?2str('b’); % allows b to be entered as a numerical value, but 


% declares b a string variable to used for 
% comparisons in logic statements. 


% Used to set the range of possible dofs for the user to choose when making changes 
% to his structure. 

% 

disp(‘Is there already a reduced FRF matrix in the correct format available,') 

disp(or does the presynthesized FRF matrix need to be generated using the dof ') 
disp(‘locations where you desire to change the structure? ') 

FRF=input('l=reduced FRF already exists 2=generate reduced FRF '); 


% Protects against user making an error in choosing how FRF is obtained. 
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%o 
if FRF~=[1 2] 
while FRF~=[1 2] 
disp(Error in choosing how FRF is obtained. Choose 1 or 2.) 
FRF=input('1=reduced FRF already exists 2=generate reduced FRF '); 
end 
end 


% Loading the presynthesized FRFs, corresponing frequency vector, and vector of dofs in 
% the eset if reduced FRF already exists. Loading of M, K, & C matrices if FRF needs to 
% generated 
%o 
if FRF== 

disp(‘Select the file which contains your presynthesized FRFs, a corresponding 
frequency ') 

disp(‘vector, and the vector of all the dofs that will be in your eset vector.’) 

pause(2) 

[hee_dofs_omega,p]=uigetfile(’*.mat',"Load [FRF], {dofs}, {omega}’); 

load (hee_dofs_omega) 


else 
disp(Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
[M_K_C,p]=uigetfile(’*.mat’,"Load [M], [K], [C]'); 
load (M_K_C) 
dofs=1:ndof; 
end 


excite=input(‘How is the structure excited? l=system 2=base '); 
% Protects against user making an error in choosing the type of excitation. 


© 
if excite~=[1 2] 
while excite~=[1 2| 
disp(‘Error in choosing type of excitation. Choose 1 or 2.’) 
excite=input(‘How is the structure excited? 1=system 2=base '); 
end 
end 


feexcite——2 

% Designation where base excitation is located. 

% 

bset=input(‘At what dof(s) are your structure excited? ‘); 
% bset=1; 


% Protects against user making an error in choosing the location of the base excitation. 
Oo 
for b_dof=1:length(bset) 
if bset(b_dof)~=dofs 


while bset(b_dof)~=dofs 
fprintf(Error in choosing location of base excitation at dof %g.\n’.... 


143 


bset(b_dof)) 
disp(‘Either the dof chosen does not exist on your structure, or 1s not ') 
disp(‘included in your list of dofs for the eset vector. *) 
bset(b_dof)=input(fprintf(‘Choose again the dof for base exitation %g “,... 
b_dof)); 
end 
end 


end 


% Designation of the spring constant which connects the substructure to the base which 
is Moving. 
% and protection against user making an error in choosing the value of the constant(s). 
Jo 
for b_spring=1:length(bset) 
kb0(b_spring)=input(fprintf(‘input the value of the spring constant which connects dof 
Jog to the base ',bset(b_spring))); 
kb0_zero=find(kb0==0); 
if length(kb0)~=b_spring | kbO_zero~=[] 
while length(kb0)~=b_spring | kbO_zero~=[] 
fprintf{(‘You made an error in entering the value for dof %g base 
exCitation spring constant.\n',bset(b_spring)) 
kb0(b_spring)=0; 
kbO(b_spring)=input(fprintf(‘Re-input the value of the spring constant 
which connects dof %g to the base ',bset(b_spring))); 
kb0O_zero=find(kb0==0); 


end 
end 
end 
cb=zeros(1,length(bset)); % initializes the damper constant which connects the 
% substructure to the base which is moving, to zero. 
end 


changek=input(‘Would you like to make suffness modifications to the structure? 1=yes 
2=no ‘); 
while changek==1 


% Designation of the spring change and protection against user making an error in 
choosing the 

% value of the constant(s). 

% 
chk=chk+1; % counter for determining # changes to stiffness matrix. 
Ik(chk)=input(‘input value of added spring (Ibs/in) ‘); 
if length(lk)~=chk 

while Jength(k)~=chk 
disp(‘You made an error entering the value for the stiffness change.’) 
Ik(chk)=0; 
Fi Ik(chk)=input(‘Re-input value of added spring (Ibs/in) ‘); 
en 
end 
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% User selection of where stiffness changes occur. Protects against user making an 
Zoerror in choosing 
% the locations of spring changes. 
%o 
kdof1(chk)=input(‘input constrained dof where 1st end of added spring is applied '); 
if kdof1(chk)~=dofs 
while kdof1(chk)~=dofs 
disp(Error in choosing location of 1st end of added spring. Either the dof 
chosen does’) 
disp(‘not exist on your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ' ) 
kdofi(chk)=input(‘Re-input constrained dof where 1st end of added spring is 
applied ’); 
end 
end 


kdof2(chk)=input(‘input dof where 2nd end of added spring is applied, b for base, or 0 
for ground ‘'); 
if kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~='b' 
while kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~="'b' 
disp(‘Error in choosing location of 2nd end of added spring. Either the dof 
chosen does’) 
disp(‘not exist on your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ‘ ) 
kdof2(chk)=input('Re-input dof where 2nd end of added spring is applied, b for 
base, or 0 for ground '); 
end 
end 


if excite==1 & kdof2(chk)=="'b' 
while excite==1 & kdof2(chk)=="b' 
disp(Error in choosing location of 2nd end of added spring. You did not 
designate this’) 
disp(‘as a base excitation structure. Therefore b is not a dof option.’) 
kdof2(chk)=input(‘Re-input dof where 2nd end of added spring is applied, b for 
base, or 0 for ground '); 
end 
end 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
7% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should be handled when asked where the structure is excited. 
%o 
if excite==2 & kdof1(chk)~=bset & kdof2(chk)=="b' 
while kdof1(chk)~=bset & kdof2(chk)=="b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(‘If this is what you desire, go back and include this dof as an excitation 
point.’) 
model_change=input(‘Was this change correct? 1=yes 2=no '); 


if model_change==1 
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error('Go back and change your basic model.') 
else 
kdof1(chk)=input('Re-input constrained dof where 1st end of added spring is 
applied ’); 
kdof2(chk)=input(‘Re-input dof where 2nd end of added spring is applied, b 
for base, or 0 for ground '); 
end 
end 
end 


% Provision for base spring constant changes to not be included in the cset matrix. 


Oo 
if excite==2 & (find(kdof1(chk)==bset))~=[] & kdof2(chk)=='b' 
cset=cset; 
else 
if kdof1(chk)~=cset % comparison of change dof with the c-set 
cset=[cset kdof1(chk)]; % matrix and formation of new c-set matrix. 
end 
if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 
end 


changek=input(‘Are there any other stiffness modifications to the structure? 1=yes 2=no 
)3 
end 
changec=input(‘Would you like to make damping modifications to the structure? l=yes 


2=no '); 
while changec==1 


% Designation of the damper change and protection against user making an error in 
choosing the 
% value of the constant(s). 
Jo 
chc=chc+1; % counter for determining # changes to damping matrix. 
Ic(chc)=input(‘input value of added damper (Ibs-sec/in) '); 
if length(ic)~=chec 
while length(Ic)~=chc 
disp(‘You made an error entering the value for the damper change.’) 
Ic(chc)=0; 
lc(chc)=input(‘Re-input value of added damper (Ibs-sec/in) '); 
end 
end 


% User selection of where damper changes occur. Protects against user making an error 
% in choosing the locations of damper changes. 

Jo 

cdof1(chc)=input(‘input constrained dof where 1st end of added damper is applied ‘'); 

if cdof1(chc)~=dofs 
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while cdof1(chc)~=dofs 
disp(Error in choosing location of 1st end of added damper. Either the dof 
chosen does’) 
disp(‘not exist on your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ' ) 
cdof1(chc)=input(‘Re-input constrained dof where 1st end of added damper is 
applied '); 
end 
end 


cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b for base, or 0 
for ground '‘); 
if cdof2(chc)~=dofs & cdof2(chc)~=0 & cdof2(chc)~='b' 
while cdof2(chc)~=dofs & cdof2(chc)~=0 & cdof2(chc)~="'b' 
disp(‘Error in choosing location of 2nd end of added damper. Either the dof 
chosen does’) 
disp(‘not exist on your structure, or is not included in your list of dofs for the’) 
disp(‘eset vector. ' ) 
cdof2(chc)=input(‘Re-input dof where 2nd end of added damper is applied, b for 
base, or 0 for ground '); 
end 
end 


if excite==1 & cdof2(chc)=='b' 
while excite==1 & cdof2(chc)=='b' 
disp(Error in choosing location of 2nd end of added damper. You did not 
designate this’) 
disp(‘as a base excitation structure. Therefore b is not a dof option.) 
cdof2(chc)=input(‘Re-input dof where 2nd end of added damper is applied, b for 
base, or 0 for ground '); 
end 
end 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from these dof to base constitutes an addition 
of 
% excitation points. This should be handled when asked where the structure 1s excited. 
% 
if excite==2 & cdof1(chc)~=bset & cdof2(chc)=="b' 
while cdof1(chc)~=bset & cdof2(chc)=='b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(‘If this is what you desire, go back and include this dof as an excitation 
point.’) 
model_change=input("Was this change correct? 1=yes 2=no ‘); 


if model_change==1 
error('Go back and change your basic model.’) 
else 
cdof1(chc)=input(‘Re-input constrained dof where 1st end of added damper is 
applied '); 
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cdof2(chc)=input('Re-input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground ‘); 
end 
end 
end 


% Provision for base damper constant changes to not be included in the cset matrix. 
% 
if excite==2 & (find(cdof1(chc)==bset))~=[] & cdof2(chc)=="b' 
cset=cset; 
else 
if cdof1(chc)~=cset % comparison of change dof with the c-set 
cset=[cset cdofl(chc)]; % matrix and formation of new c-set matrix. 
end 
if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 


. changec=input(‘Are there any other damping modifications to the structure? 1=yes 2=no 
end 

chan gem=input("Would you like to make a mass modification to the structure? l=yes 2=no 
while changem=—=1 


% Designation of the mass change and protection against user making an error in 
choosing the 
% value of the constant(s). 
Jo 
chm=chm+1; % counter for determining # changes to mass matrix. 
Im(chm)=input(‘input value of added mass (Ibs-sec’2/in) '); 
if length(lm)~=chm 
while length(lm)~=chm 
disp("You made an error entering the value for the mass change.’) 
Im(chm)=0; 
Im(chm)=input(‘Re-input value of added mass (Ibs-sec’2/in) '); 
end 
end 


% User selection of where mass changes occur. Protects against user making an error in 
% choosing the locations of mass changes. 
% 
mdof(chm)=input(‘input constrained dof where mass is added '); 
if mdof(chm)~=dofs 
while mdof(chm)~=dofs 
disp(Error in choosing location of mass change. Either the dof chosen does 
not’) 
disp(‘exist on your structure, or is not included in your list of dofs for the’) 
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disp(‘eset vector. ' ) 
mdof(chm)=input(‘input constrained dof where mass is added '); 


end 

end 

if mdof(chm)~=cset % comparison of change dof with the c-set matrix 
cset=[cset mdof(chm)]; % and formation of new c-set matrix. 

end 


changem=input(‘Are there any other mass modifications to the structure? 1=yes 2=no '); 
end 


% Handling of spring added to ground to eliminate zero from cset 
% 

ground = find(cset==0); 

cset(ground) = []; 


disp(‘The following is a list of dofs were you have made changes to your structure.') 
disp(‘This represents your {cset} vector:’) 
disp(cset) 


% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 


Oo 
if FRF==1 
disp(‘Since you inputed your own FRF matrix, your iset is chosen as all dofs remaining 


disp(‘in your eset/dof vector after extracting the cset vector.') 

iset=dofs; 

for a=1:length(cset) 

extract(a)=find(cset(a)==dofs); 

end 

iset(extract)=[]; 
else 

iset=input(‘What dofs other than those where change occured, are you interested in? ‘'); 
end 


% Formation of zset (c U b) and eset (i U c U b) matrices. 
% 
zset=[cset bset]; eset=[iset cset bset]; 


[del_Kc,del_Mc,kb] = statdelta(cset,bset,kdof1,kdof2,lk,mdof,lm,kbQ); 
save statchange M K C ndof del_Kc del_Mc kb iset cset bset zset eset excite 
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statdelta.m 


% The purpose of this function is to form the modification matrices 
Oo 


function [del_Kc,del_Mc,kb] = statdelta(cset,bset,kdof1 ,kdof2,1k,mdof,lm,kb0) 
% Formation of [del_Kc], [del_Cc], & [del_Mc] and updating of kb & cb 


Jo 


del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Mc=zeros(length(cset)); % to zero and sets their sizes 
% as [cset x cset] 
matrix 
kb=kb0; % initializes the base excitation spring constant 
% to the initial value inputed when structure was 
% formed. 
if kdof1~=[] 


for a=1:length(kdof1) 
if (find(kdof1(a)==bset))~=[] & kdof2(a)=='b' 
base=find(kdof1(a)==bset); 
kb(base) = kb(base) + lk(a); 


else 
indexk 1=find(kdof1(a)==cset); 9% finds where in cset matrix 
indexk2=find(kdof2(a)==cset); 9% the change in k occurs 
del_Kc(ndexk1,indexk 1)=del_Kc(indexk1,indexk1) + lk(a); 
del_Kc(indexk1,indexk2)=del_Kc(Gndexk1,indexk2) - lk(a); 
del_Kc(indexk2,indexk1 )=del_Kc(indexk2,indexk1) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 
end 
end 
end 
if mdof~=[] 
for a=1:length(mdof) 
indexm=find(mdof(a)==cset); % finds where in cset matrix 


% the change in M occurs 
del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + lm(a); 


end 
end 
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statsynth.m 


function [z_stat,Fred,Kred,Mred,Hstar0_dp] = statsynth(hee,kb,omega,eset,iset,bset.... 
del_Kc,del_Mc) 


% This function is designed to calculate the new static displacement using the synthesis 
% method. 


% 
pull=[); 
for w=1:length(omega) 
count=0; % initializes the counter which keeps track of 
% which column of the hee matrix is to be 
% placed into the refomation of matrix [Hee] 
for col = 1:length(eset) % reforms the [Hee] lower 
for row = col:length(eset) % triangle matrix 
count=count + 1; 
Hee_lowtri(row,col)=hee(w,count); 
end 
end 


% reforms the original symmetric [Hee] matrix 

% 

Hee=Hee_lowtn; 

[sym_row,sym_col] = find(Hee_lowtri==0); 

for n=1:length(sym_row) 
Hee(sym_row(n),sym_col(n))=Hee(sym_col(n),sym_row(n)); 

end 


e=1:length(eset); 
c=length(iset)+1:length(eset); %-length(bset); 
% bvd=length(eset)-length(bset)+1:length(eset); 
% z=length(iset)+1:length(eset); 


% Rearrangement of [Hee] to form [Hec], [Hcc], [Hce] if the structure 
% is system excited. 
Jo 
% ii excite=—=1 
Hec=Hee(e,c); 
Hce=Hee(c,e); 
Hcc=Hee(c,c); 
% end 


% Formation of [del_Zc], [del_Zb] & [del_Z] 
% 
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del_Zc=del_Kc-omega(w)*2*del_Mc; 


del_Zb=diag(kb); % diag isusedintheevent >1 base excitation 


% is present 


del_Z=[ del_Zc zeros(size(del_Zc, 1),size(del_Zb,2)) 


end 


zeros(size(del_Zb, 1),size(del_Zc,2)) del_Zb]; 


% Structure excitation Frequency synthesis equation 


oO 
H_star=Hee-Hec*del_Z*inv(eye(size(Hcc, 1)) + Hcec*del_Z)*Hce; 


% Checks to see if there are any redundant dofs in the eset due to changes and 
% bset occurring at the same dofs. If there is redundancy, it is located in 

% eset, and these rows and columns are deleted from H_star. This is done 

% because the redundant dof column(s) and row(s) in H_star creates a singular 
% matrix which is not invertable and Kred cannot be determined. 


O 
for u=1:length(bset) 
locate=find(bset(u)==eset); 
if length(locate) > 1 
pull=[pull locate(2)); 


end 

end 
H_star(pull,:)=(]; 
H_star(:,pull)=[]; 
if w==1 

Kred=inv(H_star); 

HO=H_ star; 
end 
if w==2 

H1=H_ star; 
end 
if w==3 

H2=H star; 
end 
if w==4 

H3=H_ star; 
end 


Hstar0_dp=(-H3+4*H2-5*H1+2*H0)./(omega(2)-omega(1))“2; 
Mred=.5*Kred*Hstar0_dp*Kred; 


ga=ones(size(Mred,1),1)*(-386.4); 
Fred=Mred* ga; 
z_stat=Kred\Fred; 
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APPENDIX C. TIME DOMAIN SYNTHESIS COMPUTER CODES 


The following is a list and a brief description of the main MATLAB 
computer codes that were written in order to perform the time domain 
synthesis calculations. 

e tsynconv.m - performs the time domain synthesis on a structure which 
experiences a base excitation. The dynamic responses are compared to the 
responses calculated using the classical convolution integral method. 

e tsynode45.m - performs the time domain synthesis on a structure which 
experiences a base excitation. The dynamic responses are compared to the 
responses calculated using the classical direct integration method. 

e tsynmain.m - the main program which calls the modularized time 
synthesis programs and does the post processing of the results. 

e change.m - loads the baseline structure to be modified, allows for 
modifications to the baseline structure, and calls the function delta.m 

e delta.m - forms the modification matrices which will be used to form 
impedance matrices or determine coupling forces. 

e modal.m - solves for the structure’s natural frequencies, mode shapes and 
other modal matrices and vectors. 

e impmodal.m - calculates the structure’s IRFs. 


e buildA.m - builds the quadrature matrices. 
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e fBlastForcing.m - inputs the base excitation as a function of time. 

e timesynth.m - performs the synthesis on the baseline structure and 
returns the synthesized transient responses. 

e time sift.m - selects the dynamic response that the user wants to perform 
post processing on. 


The full codes are contained on subsequent pages. Any codes that were 
previously presented will not be repeated. 
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tsynconv.m 


% The purpose of this program is to perform the time domain synthesis on a structure and 
% compare the dynamic responses to those calculated using a convolution integral method. 
Jo 


clear 


disp(‘Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
[M_K_C,p]=uigetfile(*.mat',"Load [M], [K], [C]'); 
load (M_K_C) 
dofs=1:ndof; 
C=00*C; 


[wn,phi,zeta,Mmodal,phi_norm] = modal(M,K,C); 


changek=2; changec=2; changem=2; % initializes logic variables 
% to ‘no’ status 


flag k=0; flag c=0; flag_m=0; 


chk=0; chc=0; chm=0; % initializes counter for # of changes 
% to matrices 


cset=0; % initializes cset matrix to zero to avoid null index error in 
% handling the spring to ground if no changes are made to 
% the structure. 


Hey_star=[]; Pvc]; % (re)initializes these matrices to an empty matrix 
% to avoid matrix size differences each time the 
% model is changed. 


meme —K;  Mic=M: CC c=C; % initializes change matrices which will be used in 
% the classical method, to the original matrices. 


Kb=zeros(size(K,1)); % initializes the spring, and damper to base matrices 
Cb=zeros(size(C,1)); % to zero. 


K_vec = zeros(ndof, 1); 
b=num2str(’'b’); % allows b to be entered as a numerical value, but 
% declares b a string variable to used for 


% comparisons in logic statements. 
clear MK C 


% Designation where base excitation is located, and the spring constant which 
% connects the substructure to the base which is moving. 
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% 

bset=input('At what dof(s) are your structure excited? ‘); 

for b_spring=1:length(bset); 

kb(b_spring)=input(fprintf(‘input the value of the spring constant which connects dof 

Jog to the base ',bset(b_spring))); 

end 

cb=zeros(1,length(bset)); % initializes the damper constant which connects the 
% substructure to the base which is moving, to zero. 


% Correcting the classical [K] & [Kb] matrices to include the spring element 
% connecting the substructure to the base(s). 

% 

for b_spring=1:length(bset); 


K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+kb(b_spring); 
Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring), bset(b_spring)) + 
kb(b_spring); 
K_vec(bset(b_spring))=K_vec(bset(b_spring)) + kb(b_spring); 
end 


changek=input("'Would you like to change your stiffness matrix? 1=yes 2=no '); 
while changek==1 
chk=chk+1; % counter for determining # changes 
flag k=1; 
lk(chk)=input(‘input value of added spring (Ibs/in) '); 
kdof1(chk)=input(‘input constrained dof where 1st end of added spring is applied ‘); 
kdof2(chk)=input(‘input dof where 2nd end of added spring is applied, b for base, or 
0 for ground '); 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 


O 
if kdof1(chk)~=bset & kdof2(chk)=="b' 
while kdof1(chk)~=bset & kdof2(chk)=="b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(‘If this is what you desire, go back and include this dof as an excitation 
point.') 
model_change=input(‘Was this change correct? l=yes 2=no '); 


if model_change== 
error('Go back and change your basic model.’) 
else 
kdof1(chk)=input(‘input constrained dof where 1st end of added 
damper is applied '); 
kdof2(chk)=input(‘input dof where 2nd end of added damper is applied, b 
for base, or O for ground '); 
end 
end 
end 
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% updates the base spring constant and does not count these changes 
% in the cset matrix. 


O 
if (find(kdof1(chk)==bset))~=[] & kdof2(chk)=='b' 
base=find(kdof1(chk)==bset); 
kb(base) = kb(base) + lk(chk); 
K_vec(kdof1(chk))=K_vec(kdof1(chk)) + Ik(chk); 


else 
if kdof1(chk)~=cset % comparison of change dof with the c-set 

cset=[cset kdof1(chk)]; % matrix and formation of new c-set matrix. 
en 
if kdof2(chk)~=cset 

cset=[cset kdof2(chk)]; 
end 

end 


% Classical method of reformulatting the [K] matrix. 


K_c(kdof1(chk),kdof1(chk))=K_c(kdof1(chk),kdof1(chk)) + Ik(chk); 
if kdof2(chk)~=0 & kdof2(chk)~="b' 
K_c(kdof1 (chk), kdof2(chk))=K_c(kdof1(chk),kdof2(chk)) - 1k(chk); 
K_c(kdof2(chk),kdof1(chk))=K_c(kdof2(chk),kdof1(chk)) - Ik(chk); 
K_c(kdof2(chk),kdof2(chk))=K_c(kdof2(chk),kdof2(chk)) + lk(chk); 
end 


% Classical method of reformulatting [Kb] matrix. 


O 
if kdof2(chk)=='b' 
kbb=Ik(chk); 


Kb(kdof1(chk),kdof1(chk))=Kb(kdof1 (chk), kdof1(chk)) + kbb; 
end 


changek=input(‘Are there any other changes to the stiffness matrix? l=yes 2=no ‘); 


end 


changec=input("Would you like to change your damping matrix? l=yes 2=no ‘); 
while changec==1 
che=chc+1; % counter for determining # columns in mapping matrix. 
flag_c=1; % mapping matrix flag 
Ic(chc)=input(‘input value of added damper (Ibs-sec/in) ‘); 


cdof1(chc)=input(‘input constrained dof where 1st end of added damper is applied ‘); 


cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b for base, or 
O for ground _'); 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 


‘O 
if cdof1(chc)~=bset & cdof2(chc)=="b' 
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while cdof1(chc)~=bset & cdof2(chc)=='b' 
disp(‘Connecting this dof to the base is changing the basic model.') 
disp(‘If this is what you desire, go back and include this dof as an excitation 


point.’) 

model_change=input('Was this change correct? l=yes 2=no '); 

if model_change==1 
error('Go back and change your basic model.’) 

else 
cdof1(chc)=input(‘input constrained dof where 1st end of added damper is 

applied ‘); 
cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b 
for base, or O for ground '); 
end 
end 
end 


% Updates the base damper constant and does not count these changes 
% in the cset matrix. 
%o 
if (find(cdof1(chc)==bset))~=[] & cdof2(chc)=="b' 
base=find(cdof1(chc)==bset); 
cb(base) = cb(base) + Ic(chc); 
C_vec(cdof1(chc))=C_vec(cdofl (chc)) + Ic(chc); 
else 
if cdof1(chc)~=cset % comparison of change dof with the c-set 
cset=[cset cdofl(chc)]; | % matrix and formation of new c-set matrix. 
end 
if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 


% Classical method of reformulatting the [C] matrix. 


O 
C_c(cdof1(chc),cdof1(chc))=C_c(cdof1(chc),cdof1(chc)) + Ic(chc); 
if cdof2(chc)~=0 & cdof2(chc)~="b' 
C_c(cdofl(chc),cdof2(chc))=C_c(cdof1(chc),cdof2(chc)) - Ic(chc); 
C_c(cdof2(chc),cdof1(chc))=C_c(cdof2(chc),cdof1(chc)) - le(chc); 
C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + Ic(chc); 
end 


% Classical method of reformulating [Cb] matrix. 
% 
if cdof2(chc)=="b' 

cbb=Ic(chc); - 

Cb(cdof1(chc),cdof1(chc))=Cb(cdof1 (chc),cdof1(chc)) + cbb; 
end 


changec=input(‘Are there any other changes to the damping matrix? l=yes 2=no ‘); 
end 
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changem=input(‘Would you like to change your mass matrix? 1=yes 2=no '); 
while changem==1 
chm=chm-+1; 
flag_m=1; 
Im(chm)=1nput(‘input value of added mass (lbs-sec42/in) '); 
mdof(chm)=input(‘input constrained dof where mass is added '); 


% comparison of change dof with the c-set matrix and formation of 
J new c-set matrix. 


O 
if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 


% Classical method of reformulatting the [M] matrix. 
% 
M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + lm(chm); 


changem=input(‘Are there any other changes to the mass matrix? 1=yes 2=no '); 
end 


% Handling of spring added to ground 
% 


Oo 
ground = find(cset==0); 
cset(ground) = []; 


% Formulation of iset matrix which is the set of dofs other than 

% those in the cset matrix where response information is desired. 

Jo 

iset=input("What dofs other than those where change occured, are you interested in? '); 


% Formation of zset (c U b) and eset (1 U c U b) matrices. 
Yo 
zset=[cset bset]; eset=[iset cset bset]; 


% Choosing the FRF interested in 

Jo 
resp_dof=input(‘What response dof are you interested in? '); 
frf_index=find(eset==resp_dof); 


% Formation of [del_Kc], [del_Cc], [del_Mc], [del_Zc], [del_Zb] 


% & [del_Z} 
J 
del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Cc=zeros(length(cset)); % to zero and sets their sizes 
del_Mc=zeros(length(cset)); % as [cset x cset] matrix 
if flag_k== 
for a=1:chk 
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if (find(kdof1(a)==bset))~=[] & kdof2(a)=="b' 
indexk1=[]; % provides for base changes to not be 
indexk2=[J; % included in del_Kc matrix 
else 
indexk1=find(kdof1(a)==cset); % finds where in cset matrix 
indexk2=find(kdof2(a)==cset); 9% the change in k occurs 
end 


del_Kc(indexk1,indexk1)=del_Kc(indexk 1,indexk1) + lk(a); 
del_ Kc(indexk1,indexk2)=del_Kc(indexk]1,indexk2) - lk(a); 
del_ Kc(indexk2,indexk1)=del_Kc(indexk2,indexk1) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 


end 
end 
if flag_c==1 
for a=1:chce 
if (find(cdof1(chc)==bset))~=[] & cdof2(a)== 
indexc 1=[]; 
indexc2={[]; 
else 


indexc1=find(cdofl(a)==cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); 9% the change in C occurs 
end 


del_Cc(indexcl,indexc 1)=del_Cc(indexc1,indexc1) + Ic(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexc1,indexc2) - Ic(a); 
del_Cc(indexc2,indexc 1 )=del_Cc(indexc2,indexc1) - Ic(a); 
del_Cc(indexc2,indexc2)=del_Cc(indexc2,indexc2) + Ic(a); 


end 

end 

if flag_m==1 
for a=l1:chm 


indexm=find(mdof(a)==cset); 9% finds where in cset matrix 
% the change in k occurs 


del_ Mc(indexm,indexm)=del_Mc(indexm,indexm) + Im(a); 








end 
end 
del_Kb=diag(kb); % diag is used in the event >1 base 
excitation 
del_Cb=diag(cb); Jo is present 
Jo 
Jo Time Step: 
Gas 
Start_t = 0.0; 
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%del_t = 0.0005; % plate times 


%end_t =.04; 

del_t =0.0003; % beam times 
end_t =.03; 

%del_t =0.001; % spring times 
%end_t =.2; 

time = [start_t:del_t:end_t]; % Time points 

nstep = length(time); % No. Time points 

Jo 

% Calculate Impulse Response Functions : 

Jo 

syn_t0 = clock; 

[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear wn phi Mmodal zeta 

% 


% Setup and Solve Integral Equation for x2*(t): 


% 





A = zeros(nstep); 
globalA = zeros(length(zset)*size(A,1)); 


count=0; 
col = 1+count:nstep+count; 
for acnt = 1:size(himp,1); 


for icnt_rows = 2 : nstep; 
[wts] = fTrapzWts(icnt_rows); 

A(icnt_rows, l:icnt_rows) = del_t * wts .* fliplr(himp(acnt, 1 :icnt_rows)); 
end 


row = col(1)+count:col(length(col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count = count +nstep; 
if row(length(row)) == size(globalA) 
col = col + nstep; 
count = 0; 
end 


end 
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clear A 


% 


% Checks to see if there are any redundant dofs in the eset due to changes and 
% bset occurring at the same dofs. If there is redundancy, it is located in 
% eset, and these rows and columns are deleted from H_star. This is done 
% because the redundant dof column(s) and row(s) in H_star creates a singular 
% matrix which is not invertable and Kred cannot be determined. 
% 
for u=1:length(bset) 
locate=find(bset(u)==zset); 
if length(locate) > 1 
pull=[pull locate(2)]; 
redundant=[redundant locate(1)]; 
end 
end 


pull_start=pull*nstep-(nstep-1); 
pull_end=pull*nstep; 
delete_dof={]; 
for v=1:length(pull) 
delete_dof=[delete_dof pull_start(v):pull_end(v)]; 
end 
globalA (delete_dof,:)={]; 
globalA(:,delete_dof)=[]; 


ff = ones(size(globalA,1),1); 


Yo = 1.0; % Amplitude of base motion 
[Y,Ydot] = fBlastForcing(Yo,time’, 'blst’, 0); 

tol = le-4; 

dif = 100; 


for icnt = 1:300 


X = -globalA*ff; 
if icnt >= 2 
[11,jj] = max(abs(xlast1 - X)); 


dif = max(abs(xlast1(jj) - X(j))); 
end 


xlastl = X; 
if dif < tol 

disp(‘Breaking'); break 
end 


[ff] = Synth_force(X, Y, Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t,bset,cset); 


rez 


end 


syn_tl =clock; 
time_syn_time=etime(syn_t1,syn_t0) 


% 


% Classical Method to Solve for Response: 
% ns 





clas_t0 = clock; 

Xchk = zeros(length(zset),length(time)); 
[wn_c,phi_c,zeta_c,Mmodal_c,phi_norm_c] = modal(M_c,K_c,C_c); 
%zeta_c=0; 


for icnt_modes = 1 : ndof; 
zeta_mode=zeta_c(icnt_modes); 
[{mode_irf] = fModal[RF(wn_c(icnt_modes), zeta_mode, time); 
fmodal = phi_norm_c(:,icnt_modes)' * K_vec * Y; 
Xchk = Xchk + phi_norm_c(zset,icnt_modes) * fConvTrap(mode_irf,fmodal',del_t); 
end 
clas_tl =clock; 
clas_syn_time=etime(clas_tl,clas_t0) 


% Plotting of Transient Response using Classical & Frequency Synthesis Methods 


% 


resp_index=find(zset==resp_dof); 


plot(time,Xchk(resp_index(1),:),'r--',time,X((1:nstep)+(resp_index(1)-1)*nstep),'b') 


grid 

title({"Transient Time Response at dof ',int2str(resp_dof)]) 
ylabel(‘displacement (in)’) 

xlabel(‘time (sec)') 

legend(‘Classical','Synthesis') 
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tsynode45.m 


% The purpose of this program is to perform the tme domain synthesis on a structure and 
% compare the dynamic responses to those calculated using a direct integration method. 
%o 


clear 
disp(‘Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
[M_K_C,p]=uigetfile('*.mat',"Load [M], [K], [C]'); 
load (M_K_C) 
dofs=1 :ndof; 
C=.00*E: 
[wn,phi,zeta,Mmodal,phi_norm] = modal(M,K,Q); 


changek=2; changec=2; changem=2; ~% initializes logic variables 
% to ‘no’ status 


flag _k=0; flag_c=0; flag_m=0; 


chk=0; chc=0; chm=0; % initializes counter for # of changes to matrices 


cset=0; % initializes cset matrix to zero to avoid null index error in 
% handling the spring to ground if no changes are made to the 
% structure. 
Hey_star=[]; HV_c=[]; % (re)initializes these matrices to an empty matrix 
% to avoid matrix size differences each time the 
% model is changed. 


K c=K,  Mic=M- 7G ac-@: % initializes change matrices which will be used in 
% the classical method, to the original matrices. 


K_vec = zeros(2*size(K_c,1),1); 
C_vec = zeros(2*size(C_c,1),1); 


Kb=zeros(size(K,1)); % initializes the spring, and damper to base matrices 
Cb=zeros(size(C,1)); % to zero. 
=num2str(‘b’); % allows b to be entered as a numerical value, but 


% declares b a string variable to used for 
% comparisons in logic statements. 
clear MK C 


% Designation where base excitation is located, and the spring constant which 
% connects the substructure to the base which is moving. 
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Jo 
bset=input(‘At what dof(s) are your structure excited? '); 
for b_spring=1:length(bset); 
kb(b_spring)=input(fprintf(‘input the value of the spring constant which connects dof 
%og to the base ',bset(b_spring))); 
end . 
cb=zeros(1,length(bset)); % initializes the damper constant which connects the 
% substructure to the base which is moving, to zero. 


% Correcting the classical [K] & [Kb] matrices to include the spring element 
% connecting the substructure to the base(s). 


% 
for b_spring=1:length(bset); 


K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+kb(b_spring); 

Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring)) + kb(b_spring); 

K_vec(size(K_c, 1)+bset(b_spring))=K_vec(size(K_c,1)+bset(b_spring)) + kb(b_spring); 
end 


changek=input(‘Would you like to change your stiffness matrix? l=yes 2=no '); 

while changek== 
chk=chk+1; % counter for determining # columns 1n mapping matrix. 
flag_k=1; % mapping matrix flag 
lk(chk)=input(‘input value of added spring (Ibs/in) '); 
kdof1(chk)=input(‘input constrained dof where 1st end of added spring is applied ‘); 
kdof2(chk)=input(‘input dof where 2nd end of added spring is applied, b for base, or 

0 for ground '); 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 


0 
if kdof1(chk)~=bset & kdof2(chk)=='b' 
while kdof1(chk)~=bset & kdof2(chk)=='b' 
disp(‘Connecting this dof to the base is changing the basic model.’) 
disp(‘Tf this is what you desire, go back and include this dof as an excitation 
point.") 
model_change=input(‘Was this change correct? l=yes 2=no '); 


if model_change—=1 
error('Go back and change your basic model.’) 

else 

kdof1(chk)=input(‘input constrained dof where Ist end of added damper is 

applied ’); 
kdof2(chk)=input(‘input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground '); 
end 
end 
end 


Jo updates the base spring constant and does not count these changes 
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a in the cset matrix. 


if (find(kdof1(chk)==bset))~=[] & kdof2(chk)== 

base=find(kdof1(chk)==bset); 

kb(base) = kb(base) + lk(chk); 
K_vec(size(K_c, 1)+kdof1(chk))=K_vec(size(K_c,1)+kdof1(chk)) + 1k(chk); 
else 


if kdof1(chk)~=cset % comparison of change dof with the c-set 
cset=[cset kdof1(chk)]; % matrix and formation of new c-set matrix. 
end 


if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 
end 


- Classical method of reformulatting the [K] matrix. 


K _c(kdof 1(chk),kdof1(chk))=K_ aa 1(chk),kdof1(chk)) + lk(chk); 

if kdof2(chk)~=0 & kdof2(chk)~='b 
K_c(kdof1(chk),kdof2(chk))=K_c(kdof1 (chk),kdof2(chk)) - 1k(chk); 
K_c(kdof2(chk),kdof1(chk))=K_c(kdof2(chk),kdof1(chk)) - 1k(chk); 
K_c(kdof2(chk),kdof2(chk))=K_c(kdof2(chk),kdof2(chk)) + 1k(chk); 

end 


% Classical method of reformulatting [Kb] matnx. 


Oo 
if kdof2(chk)=='b' 
kbb=Ik(chk); 


Kb(kdof1 (chk), kdof1(chk))=Kb(kdof] (chk), kdof1(chk)) + kbb; 
end 


changek=input(‘Are there any other changes to the stiffness matrix? l=yes 2=no '); 


end 


changec=input("Would you like to change your damping matrix? l=yes 2=no ‘); 
while changec==1 
chc=chc+1; % counter for determining # columns in mapping matrix. 
flag_c=1; % mapping matrix flag 
Ic(chc)=input(‘input value of added damper (Ibs-sec/in) '); 


cdof1(chc)=input(‘input constrained dof where Ist end of added damper is applied ‘); 
cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b for base, or 
0 for ground '); 


% Let's the user know that he is not allowed to make dof(s) to the base changes which 
Yo are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 


oO 
if cdof1(chc)~=bset & cdof2(chc)=="b' 
while cdof1(chc)~=bset & cdof2(chc)== 
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disp(‘Connecting this dof to the base is changing the basic model.') 

disp(‘If this is what you desire, go back and include this dof as an excitation 
point.’) 

model_change=input('Was this change correct? l=yes 2=no '); 


if model_change==1 
error('Go back and change your basic model.’) 


Ee 
cdof1(chc)=input(‘input constrained dof where 1st end of added damper is 
applied ’); 
cdof2(chc)=input(‘input dof where 2nd end of added damper is applied, b 
for base, or O for ground '); 
end 
end 
end 


% Updates the base damper constant and does not count these changes 
% in the cset matrix. 
% 
if (find(cdof1(chc)==bset))~=[] & cdof2(chc)=="b' 
base=find(cdof1(chc)==bset); 
cb(base) = cb(base) + Ic(chc); 
C_vec(size(C_c,1)+cdof1(chc))=C_vec(size(C_c,1)+cdofl(chc)) + lc(chc); 
else 
if cdof1(chc)~=cset % comparison of change dof with the c-set 


cset=[cset cdofl(chc)]; | % matrix and formation of new c-set matrix. 
end 


if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 


: Classical method of reformulatting the [C] matrix. 


Cc _c(cdof1(chc),cdof1(chc))=C_ c(cdof 1(chc),cdof1 (chc)) + Ic(chc); 

if cdof2(chc)~=0 & cdof2(chc)~='b' 
C_c(cdof1(chc),cdof2(chc))=C_c(cdof (che),cdof2(chc)) - lc(chc); 
C_c(cdof2(chc),cdof1(chc))=C_c(cdof2(chc),cdof1(chc)) - Ic(chc); 


C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc); 
end 


% Classical method of reformulating [Cb] matrix. 


‘0 
if cdof2(chc)== 
cbb=Ic(chc); 


Cb(cdof1(chc),cdof1(chc))=Cb(cdof1 (chc),cdof1(chc)) + cbb; 
end 


changec=input(‘Are there any other changes to the damping matrix? 1=yes 2=no _); 
end 
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changem=input(‘Would you like to change your mass matrix? l=yes 2=no '); 
while changem==1 
chm=chm+1; 
flag_m=l1; 
Im(chm)=input(‘input value of added mass (Ibs-sec42/in) '); 
mdof(chm)=input(‘input constrained dof where mass is added '); 


% comparison of change dof with the c-set matrix and formation of 
J new c-set matrix. 


oO 
if mdof(chm)~=cset 
cset=[cset mdof(chm)}; 
end 


% Classical method of reformulatting the [M] matrix. 
% 
M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + Im(chm); 


changem=input(‘Are there any other changes to the mass matrix? 1=yes 2=no '); 
end 


% Handling of spring added to ground 
qe 


O 
ground = find(cset==0); 
cset(ground) = []; 


% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 


% 
iset=input('What dofs other than those where change occured, are you interested in? ‘); 


% Formation of zset (c U b) and eset a U c U b) matrices. 
% 
zset=[cset bset]; eset=[iset cset bset]; 


% Choosing the FRF interested in 

% 
resp_dof=input(‘What response dof are you interested in? ‘); 
frf{_index=find(eset==resp_dof); 


% Formation of [del_Kc], [del_Cc], [del_Mc], [del_Zc], [del_Zb] 


% & [del_Z] 

70 | 
del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Cc=zeros(length(cset)); % to zero and sets their sizes 
del_Mc=zeros(length(cset)); % as [cset x cset] matrix 

if flag_k== 

for a=1:chk 


if (find(kdof1(a)==bset))~=[] & kdof2(a)=='"b' 
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indexk1=[]; % provides for base changes to not be 
indexk2=[]; % included in del_Kc matrix 
else 
indexk 1=find(kdof1(a)==cset); % finds where in cset matrix 
indexk2=find(kdof2(a)==cset); 9% the change in k occurs 
end 


del_Kc(indexk1,indexk1)=del_Kc(indexk1,indexk1) + Ik(a); 
del_Kc(indexk1,indexk2)=del_Kc(indexk1,indexk2) - lk(a); 
del_Kc(indexk2,indexk1)=del_Kc(indexk2,indexk1) - Ik(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 


end 
end 
if flag_c==1 
for a=l:chc 
if (find(cdof1(chc)==bset))~=[] & cdof2(a)=="b' 
indexc1=[]; 
indexc2=[]; 
else 


indexcl1=find(cdof1(a)==cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); % the change in C occurs 
end 


del_Cc(indexc1,indexc1)=del_Cc(indexc1,indexc1) + Ic(a); 
del_Cc(indexc1,indexc2)=del_Cc(indexc1,indexc2) - Ic(a); 
del_Cc(indexc2,indexc1)=del_Cc(indexc2,indexc1) - Ic(a); 
del_Cc(indexc2,indexc2)=del_Cc(indexc2,indexc2) + Ic(a); 
end 
end 


if flag_m== 
for a=1:chm 
indexm=find(mdof(a)==cset); 9% finds where in cset matrix 
% the change in k occurs 


del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + Im(a); 





end 
end 

del_Kb=diag(kb); % diag is used in the event >1 base excitation 
del_Cb=diag(cb); % is present 

% 

% Time Step: 

Oa: 

Start_t = 0.0; 

%del_t =0.00025: % plate times 

%end_t =.02; 
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del t =0.0003; % beam times 


end_t =.03; 

%del_t =0.001; % spring times 
Fendt =.2; 

time = [start_t:del_t:end_t]; % Time points 

nstep = length(time); % No. Time points 

%o 

% Calculate Impulse Response Functions : 

ly minis esse eta ee 

syn_t0O = clock; 

[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear wn phi Mmodal zeta 

Jo 


% Setup and Solve Integral Equation for x2*(t): 


% 





A = zeros(nstep); 
globalA = zeros(length(zset)*size(A, 1)); 


count=0; 
col = 1+count:nstep+count; 
for acnt = 1:size(himp, 1); 


for icnt_rows = 2 : nstep; 
[wts] = fTrapzWts(icnt_ rows); 

A(icnt_rows, l:icnt_rows) = del_t * wts .* fliplr(himp(acnt,1:icnt_rows)); 
end 


row = col(1)+count:col(length(col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count = count +nstep; 


if row(length(row)) == size(globalA) 
col = col + nstep; 
count = OQ; 
end 
end 
clear A 
Zodisp(sprintf(‘Norm(globalA) = %5.3f ,norm(globalA))) 
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% 
ff = ones(length(zset)*nstep, 1); 


Yo = 1.0; % Amplitude of base motion 
[Y,Ydot] = fBlastForcing(Yo,time’, 'blst', 0); 

tol = le-4; 

dif = 100; 


for icnt = 1:300 
X = -globalA*ff; 
if icnt >= 2 
[ii,])] = max(abs(xlast1 - X)); 


dif = max(abs(xlast1(jj) - XQj))); 
end 


xlastl = X; 
if dif < tol 


disp(‘Breaking'); —_ break 
end 


[ff] = Synth_force(X, Y, Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t, bset,cset); 
end 


syn_tl =clock; 
time_syn_time=etime(syn_t1,syn_t0) 


% 
% Classical Method to Solve for Response: 


(@) GP AAG I ALD fA OG IOGEAR EG OLE LLG CLOG LLG ALE ALL GEA LLG (ALGAE ALG ALG ALLE ALL LLG LEG AI IGE ARG A LEIDER 


clas_tO =clock; 

Save structure K_c C_c M_c K_vec C_vec Yo 
XchkO0=zeros(2*ndof,1); 
[tchk,Xchk]=ode45(‘scnd_to_frst’,start_t,end_t,XchkQ); 
clas_tl =clock; 

clas_syn_time=etime(clas_tl,clas_tO) 


% Plotting of Transient Response using Classical & Frequency Synthesis Methods 


% 


resp_index=find(zset==resp_dof); 


V7 


plot(tchk,Xchk(:,resp_dof(1)),'r--',tume,X((1:nstep)+(resp_index(1)-1)*nstep),'b’) 
id 


title({["Transient Time Response at dof ',int2str(resp_dof)]) 
ylabel(‘displacement (in)’) 

xlabel(‘time (sec)') 

legend(‘Classical', Synthesis’) 
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tsynmain.m 


% This program is designed to calculate the new transient responses of a system after 
% changes in the mass, stiffness, and damping matrices have been made. 
% The new responses will be calculated using the synthesis method in the time domain. 


% 
clear 


load change 

% Designation of the frequency range and step size 
Jo 

% initial stepsize end 

% Pee 
t_input =[ 0.0 0.0004 BU a 





PD 








time = t_input(1):t_input(2):t_input(3); 
nstep = length(tme); % No. Time points 


[wn,phi,zeta,Mmodal,phi_norm] = modal(M,K,C); 
clear MK C 


syn_t0O = clock; 


[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear wn phi zeta Mmodal 


[globalA] = buildA(himp,bset,zset,nstep,t_input(2)); 


Yo = 1.0; % Amplitude of forcing function 
[Y, Ydot] = fBlastForcing(Yo,time’, 'blst'’, 0); 


[X_star,icnt] = timesynth(globalA,zset, Y, Ydot,del_Kc,del_Cc,del_Mc,kb,cb.... 
nstep,t_input(2),bset,cset); 


syn_tl = clock; 
time_syn_time=etime(syn_t1,syn_t0) 


% Plotting of Transient Response using Classical & Frequency Synthesis Methods 
E idot=inpurwhat response dof are you interested in? ‘); 

[X_star_desired] = time_sift(zset,resp_dof,X_star,nstep); 
plot(time,X_star_desired,’b’) 


grid 
title({'Synthesized Transient Time Response at dof ',int2str(resp_dof)]) 
ylabel(‘displacement (in)’) 
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xlabel(‘time (sec)') 





function [himp] = impmodal(wn,phi,zeta,Mmodal,t,imp_dof) 
% This function is designed to calculate the IRF using mode shapes and 


% summing the IRF over a number of modes. 
% 


ncol=(length(imp_dof)*(length(imp_dof)+1))/2; 


% Calculation of impulse response function 


Jo 
himp=zeros(ncol,length(t)); % Initialization to zero, and sizing 
% of the freq x lower triangle matrix. 
nmodes=size(phi); b=0; % Defining the number of modes that exist and 
% initializing b which keeps track of the 
% of modes 


for b=1:nmodes 
% Elimination of unnessecary elements and rearrangement of original 
% mode shapes {phi} to form {phi_red} and [num_red] which is now an 
% imp_dof x imp_dof matrix. 
Jo 
num_red=phi(imp_dof,b)*phi(imp_dof,b)’; 
if wn(b) > le-3 % ‘Then elastic mode 


wd = wn(b) * sqrt(1 - zeta(b)“2); 
mode_irf = exp(-zeta(b)*wn(b)*t) .* sin(wd * t) / (Mmodal(b)*wd); 


elseif wn(b) <= le-3 % Rigid body mode 
mode_irf = t(Mmodal(b); 
end 
% Changes the [num_red] matrix from an imp_dof x imp_dof matrix to 
% a 1x lower triangle vector. 


% 


H_lowtn=tnl(num_red); % pulls out the lower triangle 
% and places zeros everywhere else 


symtry = find(H_lowtn==0); % finds zero values in the 
| J lower triangle matrix 


H_lowtri(symtry) = []; % deletes all values of zero and 
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% turns the remaining elements 
% into a 1 x length(lower triangle) 
% vector 


for cnt=1:length(H_lowt) 
himp_mode(cnt,:)=H_lowtni(cnt)*mode_uirf; 
end 


himp = himp + himp_mode; 
end 
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buildA.m 


function [globalA] = buildA(himp, bset,zset,nstep,del_t) 


% The purpose of this program is to build the trapezoidal quadrature matrices to be used 
% for the solution of the VIDE 
Jo 


A = zeros(nstep); 
globalA = zeros(length(zset)*size(A, 1)); 


count=0; 
col = 1+count:nstep+count; 
for acnt = 1:size(himp, 1); 


for icnt_rows = 2 : nstep; 
[wts] = fTrapzWts(icnt_rows); 

A(icnt_rows, l:icnt_rows) = del_t * wts .* flipir(himp(acnt, i :icnt_rows)); 
end 


row = col(1)+count:col(length(col))+count, 
globalA(row,col)=A; 
globalA(col,row)=A; 
count = count +nstep; 


if row(length(row)) == size(globalA) 
col = col + nstep; 
count = Q; 
end 
end 


% Checks to see if there are any redundant dofs in the eset due to changes and 
% bset occurring at the same dofs. If there is redundancy, it is located in 
% eset, and these rows and columns are deleted from H_star. 


Oo 
for u=1:length(bset) 
locate=find(bset(u)==zset); 
if length(locate) > 1 
pull=[pull locate(2)]; 
redundant=[redundant locate(1)]; 
end 
end 


pull_start=pull*nstep-(nstep-1); 
pull_end=pull*nstep; 
delete_dof=[]; 
for v=1:length(pull) 
delete_dof=[delete_dof pull_start(v):pull_end(v)]; 
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end 
globalA (delete_dof,:)=[]; 
globalA (:,delete_dof)=[]; 
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fBlastForcing.m 


function (f_of_t,fdot] = fBlastForcing(Fo,time, type, plotit); 


% 
% Usage: [f_of_t,fdot] = fBlastForcing(Fo,time, type, plotit); 
% 
% Choices: sine blst step 
% 
%  type='step’ STRING Variable 
iow 
% 
% 
% If use ‘sine’, fdot also returned. 
% 
% This function returns a forcing function which is 
% a"blast" function. 
% 
% F(t) = Fo * ( exp(-at) - exp(-bt) ) 
% 
% where a and b are constants which shape the blast, 
% and Fo is the amplitude of the blast. 
% 
% The variable "plotit" is a switch which if set = 1 will 
% cause the f(t) to be plotted, if set to anything else 
% will not plot. 
% 
% 
% Choices: sine blst step 
% type ='step;; 
i re ans 
if type == ‘bist’; 
%  disp(' Blast forcing used...') 
a = 100.0; 
b = 300.0; 


f_of_t = Fo * ( exp(-a*time) - exp(-b*time) ); 
fdot = Fo * ( -a*exp(-a*time) + b*exp(-b*time) ); 


elseif type == 'step’; 


% 


disp(' Step forcing used...’ 
f_of_t = Fo * ones(size(time)); 
fdot = [}; 


elseif type == ‘sine’; 
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%  disp(' Sine forcing used...') 
W =5; % Hertz 
f_of_t = Fo * sin(2*pi*W*time); 
fdot = Fo * (2*pi*W)*cos(2*pi* W*time); 


end; 


if plout = 1; 
figure(gcf+1) 
if type == ‘sine’; 
plot(time,f_of_t,time,fdot);grid 
else 
plot(time,f_of_t); grid 
end 
pause 
% cif 
end 


% End function. 
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timesynth.m 


function [X_star,icnt] = timesynth(globalA,zset, Y, Ydot,del_Kc,del_Cc,del_Mc,kb,cb.,... 
nstep,del_t,bset,cset) 


% This function is designed to calculate the new dynamic response using the synthesis 
% method. 


J 

ff = ones(size(globalA, 1), 1); 
tol = le-2; 

dif = 100; 


for icnt = 1:300 
X_star = -globalA*ff; 


if icnt >= 2 

[11,]j] = max(abs(xlast1 - X_star)); 

dif = max(abs(xlast1(jj) - X_star(jj))); 
end 


if icnt>=300 
li 
ay 
dif 

end 


xlastl = X_star; 


if dif < tol 
% disp(‘Breaking’); 
break 
end 
% icnt 


[ff] = 
Synth_force(X_star, Y, Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t, bset,cset); 


end 
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function [X_star_desired] = time_sift(zset,resp_dof,X_star,nstep) 

% The purpose of this program is to sort through the synthesized transient 
% responses and return the one in which the user is interested in. 

% 

resp_index=find(zset==resp_dof); 


X_Star_desired = X_star((1:nstep)+(resp_index(1)-1)*nstep); 
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APPENDIX D. OPTIMIZATION COMPUTER CODES 


The following is a list and a brief description of the main MATLAB 


computer codes that were written in order to perform the optimization of a 


shock and vibration isolation system. 


isomain.m - the main program which loads the baseline structure to be 
optimized, and allows for the designation of the design variables. Also, if 
necessary, calls various modularized functions such as modal.m, 
frfmodal.m, impmodal.m, fBlastForcing.m, buildA.m. It then allows for 
the optimization inputs, the use of the MATLAB Optimization function 
constr.m, and the post processing. 

modal.m - solves for the structure’s natural frequencies, mode shapes and 
other modal matrices and vectors. 

frfmodal.m - calculates the structure’s FRFs. 

impmodal.m - calculates the structure’s IRFs. 

buildA.m - builds the quadrature matrices. 

fBlastForcing.m - inputs the base excitation as a function of time. 


isosub.m - the subroutine program which is called by the MATLAB 


. Optimization function constr.m. This is where the optimization iterations 


occur, the design variables change, and the system’s responses are 


calculated. This program calls the functions delta.m, frfsynth.m, 
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statdelta.m, statsynth.m and timesynth.m. It then calculates the 
normalized constraint values. 

e delta.m - forms the modification matrices which will be used to form 
impedance matrices or determine coupling forces. 

e frfsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized FRFs. 

e statdelta.m - forms the modification matrices which will be used to form 
impedance matrices. 

e statsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized static displacements. 

e timesynth.m - performs the synthesis on the baseline structure and 
returns the synthesized transient responses. 


The full codes are contained on subsequent pages. Any codes that were 
previously presented will not be repeated. 
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clear 

clear global 

diary isomain.dia 
tic 


global kdof1 kdof2 cdof1 cdof2 mdof iset cset bset zset eset kbO cb0 hee omega ... 
opt_omega excite resp_dof inp_dof resp_index stat_omega stat_hee iter ... 
del_vars del_f H_desired_0 del_H_desired globalA Y Ydot t_input nstep 


Oo 7 6 A A A Ae He Ae Ae He He He 2 ee Ae ee ee Ae Ae He Ae 2 He He He he He 2 2 2 Ae ee 2 ee he He he Ae Ae 2 2 2 i 2 2 2 2 A ee ie Ae Ae A Ae Ae ae 2 2 A 2 


Jo 

Jo Loading of Structure to be Optimized 
sens 
Jo 

disp(‘Are there already reduced FRF & IRF matrices in the correct format available,') 
disp(‘or does the presynthesized FRF & IRF matrices need to be generated using the dof ') 
disp(‘locations where you desire to change the structure? ') 

FRF_IRF=input('l=reduced FRF/IRF already exists 2=generate reduced FRF/IRF '); 





% Protects against user making an error in choosing how FRF/IRF is obtained. 


Oo 
if FRF_IRF~=[1 2] 
while FRF_IRF~=[1 2] 
disp(‘Error in choosing how FRF/IRF is obtained. Choose 1 or 2.') 
FRF_IRF=input('1=reduced FRF/IRF already exists 2=generate reduced FRF/IRF '); 
end 


end 
if FRF_IRF==1 
disp(‘Select the file which contains your presynthesized FRF/IRFs, a corresponding 
frequency ') 
disp(‘and time vectors, and the vector of all the dofs that will be in your eset vector.’) 
pause(2) 


[hee_dofs_omega,p]=uigetfile(**.mat',,Load FRF/IRF’) 
load (hee_dofs_omega) 


disp(‘Also select the file which contains your presynthesized FRFs with zero damping ’) 
disp(‘for the purpose of calculating static displacement.’) 

pause(2) 

[stat_hee,p]=uigetfile('*.mat'’,"Load Static FRF’) 

load (stat_hee) 


else 


disp(‘Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
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[M_K_C,p]=uigetfile(’*.mat'’,"Load [M], [K], [C]’); 
load (M_K_C) 
dofs=1:ndof; 
kc0=K(109,109); 
end 


fy 5 Ae 2 Ae 2 Ae 2 2 Ae Ee Ae 2 ie He He Ae Ae 2c He 2 Ae 2 2c ie ee 2 2 he he 2 2 2c Ae 2 2 2 ie 2 ie 2 2 2 Ae ie 2 2c Ae 2 2c 6 2 2 ie 2 ie 2 2 ie 2 2 ae 2 A a a 


% 
J Designation of Optimization and Synthesis DOFs 

Oly exerting 

% 

kdof1=[]; kdof2=[]; % initalizes the matrices which keep 
cdof1=[}; cdof2=[{]; % track of the dofs where changes occur 
mdof={[]}; 


b=num2str('b’); % allows b to be entered as a numerical value for the 
dof 

% choices; but declares b a string variable to 
used for 

% comparisons in logic statements. 


kdofs=input(Input dof pairs for suffness changes for optimization '); 
if kdofs==[] 

kdofl=[]; kdof2={]; 
else 
kdofl=kdofs(:,1); kdof2=kdofs(:,2); 
end 


cdofs=input(Input dof pairs for damping changes for optimization ‘); 
if cdofs==[] 
cdofl=[]; cdof2=[]; 
else 
cdofl=cdofs(:,1); cdof2=cdofs(:,2); 
end 


mdof=input(‘Input dofs for mass changes for optimization '); 


Up FAIR IRC ACR ICRICR AOR II IRC IC II CR IC 


%o 
% Formation of Partitioning and Arrangement Vectors 
Oly remeron tee ttt ci Pg ee 


% 
excite=input(‘How is the structure excited? l=system 2=base '); 


% Formation of {bset} 

% 

cset=0; bset=[]; % initializes cset matrix to zero to avoid 
% null index error if no changes are made to 
% the structure and bset to empty matrix to 
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% prevent error in the event this is not a 
% base excitation and bset does not exist. 


kbO=0; cb0=0; % initializes kb and cb to zero to prevent error in 
% the event this is not a base excitation and kb & cb does 
% not exist. 


if excite==2 
% Designation where base excitation is located. 
% 
bset=input(‘At what dof(s) are your structure excited? ‘); 


% Designation of the spring and damper constants which connects the substructure to the 
% base which is moving.and protection against user making an error in choosing the 
% value of the constant(s). 
% 
for b_spring=1:length(bset) 
kbO(b_spring)=input(fprintf(‘input the value of the spring constant which connects 
dof %g to the base ',bset(b_spring))); 
kbO_zero=find(kb0==0); 
if length(kbO)~=b_spring | kbO_zero~=[] 
while length(kbO)~=b_spring | kbO_zero~=[] 
fprintf("You made an error in entering the value for dof %g base excitation 
spring constant.\n',bset(b_spring)) 
kbO(b_spring)=0; 
kbO(b_spring)=input(fprintf(‘Re-input the value of the spring constant which 
connects dof %g to the base ',bset(b_spring))); 
kbO_zero=find(kb0==0); 


end 
end 
end 
cb0=zeros(1,length(bset)); % initializes the damper constant which connects the 
% substructure to the base which is moving, to zero. 
end 


% Formation of {cset} vector 


oO 
if kdofs~=[] 
for chk=1:length(kdof1) 


% Comparison of change dof with the c-set matrix and formation of new c-set matrix. 
% Does not count changes to the base spring constant in the cset matrix. 
% 
if excite==2 & (find(kdof1(chk)==bset))~=[] & kdof2(chk)=="b' 
cset=cset; 
else 
if kdof1(chk)~=cset 
cset=[cset kdof1(chk)]; 
end 
if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
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end 
end 
end 
end 


if cdofs~=[] 
for chc=1 :length(cdof1) 


% Comparison of change dof with the c-set matrix and formation of new c-set matrix. 
: Does not count changes to the base damper constant in the cset matrix. 


if exCcite==2 & (find(cdof1(chc)==bset))~=[] & cdof2(chc)=="b' 


cset=cset; 
else 
if cdof1(chc)~=cset 
cset=[cset cdofl(chc)]; 
end 
if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 
end 
end 
if mdof~=[] 


for chm=1:length(mdof) 


% Comparison of change dof with the c-set matrix and formation of new c-set matrix. 
% 


if mdof(chm)~=cset 
cset=[cset mdof(chm)]; % and formation of new c-set matrix. 
end 
end 
end 
ground = find(cset==0); % Handling of spring added to ground to 
eliminate 
cset(ground) = []; % zero from cset 


% Formation of {iset} vector 


O 

if FRF_IRF==1 

disp('Since you inputed your own FRF/IRF matrix, your iset is chosen as all dofs 
remaining ') 

disp(‘in your eset/dof vector after extracting the cset vector.’) 

iset=dofs; 

for a=1:length(cset) 

extract(a)=find(cset(a)==dofs); 

end 

iset(extract)=[]; 
else 
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% iset=input(‘What dofs other than those where change occured, are you interested in? 


iset=[]; 
end 
% Formation of zset (c U b) and eset (i U c U b) matrices. 


J 
zset=[cset bset]; eset=[iset cset bset]; 


Tf 2 2 A He 2 Ae He He 2 ee Ae 2 ie ee Ae Ae 2 Ae he 2 26 ie 2 Ae ie 2 2 2 he 2 2 2 ie ie 2 Ae oe ee ok 2 2 he 2 2 ie 2k he oe 2 ee 2 2c 2 oe 2 ie 2 ek 2 2k 2 


% Formulation of original [H], [Himp], for Freq & Time Synthesis and Static 
J Displacement 


% 

if FRF_IRF==2 
% Designation of the frequency and time range and step size 
% 
J initial stepsize end 


Oo Pa Fat Pt Poet Paget Pt Pt Pap Pa PRP Pat Frat Pad Pap fart Pind Pag PUD PAD PD Pt Od 


freq = [ Ol 10 1e3+.01 it: 
t_input =[ 0.0 0.001 05 }; 


omega=freq(1):freq(2):freq(3); 
time = t_input(1):t_input(2):t_input(3); 
nstep = length(tme); % No. Time points 


[wn,phi,zeta,Mmodal,phi_norm,Cmodal] = modal(M,K,C); 
[hee] = frfmodal(wn,phi,zeta,Mmodal,omega,eset); 
[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 


% Information for static displacement synthesis 
% 
stat_zeta = 0* zeta; 
Stat_omega = .1:.1:.4; 
[stat_hee] = frfmodal(wn,phi,stat_zeta,Mmodal,stat_omega,eset); 


end 
clear M K C 


Off 78 AE A He AC Ae A He He A Ae A He He A A A A A RK A RAK A KA AERA KARE KKK AK KAKA A KK AK KK KK 


Jo 

% Time Synthesis 

noon 

% Builds the matrix A which uses trapezoidal weights to integrate the [RFs 
% 


O 
[globalA] = buildA(himp,bset,zset,nstep,t_input(2)); 
clear himp 
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% Forms the base excitation vectors of displacement and velocity 
% 

Yo= 1.0; % Amplitude of forcing function 

[Y, Ydot] = fBlastForcing(Yo,time’, 'blst’, 0); 


7 Ae A i Ae Ae Ae Ae He 2 He He Ae He 2 2 ee Ae he 2 Ae he 2 ee Ae He ae 2 2c Ae ee 2 2c Ae ae 2 2c he he He ae 2 he 2 2 he he 2 2c Ae fe 2c Ae 2h he he ae 2c he 2 2c 2c 2 2 2c 2c 2k 


% 
% The Optimization 
Oy mamta acs eerie eee 


% Designates the frequency over which the optimization will be conducted 
% 
if FRF_IRF==2 
disp(Since you generated the FRF/IRFs through this optimization program, your FRF’) 
disp(‘frequencies and IRF times are automatically chosen as the freq. and time range ') 
disp(‘for your optimization.) 
opt_omega=omega; 
else 
opt_freq1=input(‘What is the initial frequency over which to evaluate this problem’); 
opt_freq2=input('What is the final frequency over which to evaluate this problem’); 
disp(‘The frequency increment is automatically chosen to match the FRF freq. 
increment’) 
opt_omega=opt_freq 1:omega(2)-omega(1):opt_freq2; 


end 


% Selecting the FRF that is to be minimized. 
% 
resp_dof=input(‘What response dof are you interested in? ‘); 


if excite==2 
inp_dof=[]; 


inp_dof=input(‘What input dof are you interested in? '); 
end 
resp_index=find(zset==resp_dof); 


% Setting the number of of optimization variables, their initial values, and 
% their upper and lower bounds. 

% | 

numoptvar=input(‘How many optimization variables do you have? '); 


delO=[1 1 1 1]; 


% del_kp del_kc del_mp 

% Fe smacks ee 

% del_kb del_kc del_cb del_cc 

% mnovansmabdet eld. SiR =, Ay 

delib= [ -4 -4 0 0 I; 
delub= [ 5 5 5 5 }; 
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%o Initializing the matirices which will keep track of how the design variables 


% and objective function changes for every optimization iteration for the 
% purpose of plotting. 

% 

del_vars =[]; del_f=[]; del_H_desired = []; 

iter=0; % Sets the number of iterations counter at zero 

% Calling of constr.m routine to perform optimization 

% on subroutine iso_sub.m 

% 


options(1)=1; 
options(14)=200; 
toc 


tic 
[del]=constr(‘isosub’,del0,options,dellb,delub); 


Of 26 2 Ae 2 Ae Ae ee Ae 2 Ae Ae 2 Ae Ae ie Ae 2 ie Ae 2 ee 2 2 Ae Ae Ae 6 2 Ae 2 2 2 Ae 2 Ae Ae Ae Ae ie Ae 2 2 Ae 2 Ae Ae ee A Ae 2 2 Ae 2 2 AG 2 A A 26 oe 2 AK A A 


% OUTPUT 


% Peet Peat Pet Pet Ped Pd Poet Peat Pet Pt 


load opt_data 


final_del_kb = del_vars(length(del_vars), 1); 
final_del_kc = del_vars(length(del_vars),2); 
final_del_cb = del_vars(length(del_vars),3); 
final_del_cc = del_vars(length(del_vars),4); 


disp("The recommended change in base isolator stiffness (Ib/in) is ---> '),final_del_kb 
disp("The recommended change in computer isolator suffness (lb/in) is ---> '),final_del_kc 
disp("The recommended change in base isolator damping (1b/in/s) is ---> '),final_del_cb 
disp("The recommended change in computer isolator damping (Ib/in/s) is ---> 
),final_del_cc . . 

disp(‘Constraints:’) 

disp(g) 


iterations=1:iter; 
% Plot of FRF to observe how it changes during the optimization 


% 

figure(1) 
semilogy(omega,abs(del_H_desired(1,:)),'--r',omega,abs(del_H_desired(2,:)),'m’) 
utle(‘(Optimized Frequency Response Function [H*]') 

ylabel(‘FRF Amplitude (unitless)’) 

xlabel(‘Frequency (rad/s)') 

legend('Before’,'After’) 


figure(2) 
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subplot(2,1,1) 
plot(iterations,del_vars(:,1),'--r',iterations,del_vars(:,2),'m’) 
title(‘Change in Variables vs. Iterations’) 
ylabel(‘Isolator Stiffness (1b/in)’) 
legend(‘del_kb’,'del_kc’) 

subplot(2,1,2) 
plot(iterations,del_vars(:,3),'--r',iterations,del_vars(:,4),'m’) 
title(‘Change in Variables vs. Iterations’) 
xlabel(‘# of iterations’) 
ylabel(‘Isolator Damping (lb-s/in)’) 
legend(‘del_cb’,'del_cc’) 


figure(3) 

subplot(2,1,1) 
plot(iterations,del_vars(:,1)+kb0(1),'--r',iterations,del_vars(:,2)+kc0,'m') 
title((Change in Isolation Constants vs. Iterations’) 
ylabel(‘Isolator Suffness (1b/in)') 
legend(‘k_base’,'k_computer’) 

subplot(2,1,2) 
plot(iterations,del_vars(:,3),'--r',1terations,del_vars(:,4),'m') 
title((Change in Isolation Constants vs. Iterations’) 
xlabel('# of iterations’) 
ylabel(‘Isolator Damping (lb-s/in)’) 
legend(‘c_base’,'c_computer’) 


figure(4) 

semilogy(iterations,del_f,'r’) 

title(‘Change in Objective Function vs. Iterations’) 
ylabel('Magnitude’) 

xlabel(‘# of iterations’) 


figure(5) 
plot(time,X_star((1:nstep)+(resp_index(1)-1)*nstep),'b’) 
id 


gri 
title([‘Synthesized Transient Time Response at dof ',int2str(resp_dof)]) 
ylabel(‘displacement (1n)’) 

xlabel(‘time (sec)’) 


figure(6) 
plot(time,Xaccel_star/386.4,'m’) 


grid 

title((‘Synthesized Transient Time Response at dof ',int2str(resp_dof)]) 
ylabel(‘acceleration (g"s)’) 

xlabel(‘time (sec)’) 


Ty BRK RK KER E KR KKERKKRKK KKKKKKKRK KKK KKK KEK AKA KAEKAKAKK KK KKKKKK KAKA KEK KKK 
toc 

diary off 

% end Iso_main.m 
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global kdof1 kdof2 cdofl cdof2 mdof iset cset bset zset eset kbO0 cbO hee omega... 
opt_omega excite resp_dof inp_dof resp_index stat_omega stat_hee iter del_vars... 
del_f H_desired_0 del_H_desired globalA Y Ydot t_input nstep 


function [f,g]=iso_sub(del) 


iter=iter + 1; % counts number of iterations 


% Initialization and assigning of the optimization variables to changes which are made to 
% the structure. Spring and damper constants are scaled to even out the domain. 


oO 
lk=0; Ic=0; Im=0; 


del_kb=1e3*del(1); del_kc=1e3*del(2); 
del_cb=1*del(3); del_cc=1*del(4); 
1k(1:4)=del_kb*ones(1,4); 1k(5)=del_kc; 
1k(6:9)=-1e3*ones(1,4); 

lc(1:4)=del_cb*ones(1,4); lc(5)=del_cc; 


Of 95 6 He AA Ae He He Ae He He He HH He He He A he 6 He He 2 eC A ie fe He he 26 2 2 2 2 2 2 2 2 he he 2 2h hee 2 2 2 2 2 2c ee 2 2 ie he 2k 2k 2 Ae Ae ee 2k 
% 

% Formulation of FRF and Calculation of mean square acceleration 

% Fa a a a a a Raa 


(del_Kc,del_Cc,del_Mc,kb,cb] = 
delta(cset, bset,kdof1 ,kdof2,lk,cdof1,cdof2,lc,mdof,lm,kb0,cb0); 


[h_star] = frfsynth(hee,kb,cb,omega,excite,eset,iset,bset,del_Kc,del_Cc,del_Mc); 


[H_desired] = frf_sift(eset,resp_dof,inp_dof,excite,h_star); 
cut_freq=[find(omega<opt_omega(1)) find(omega>opt_omega(length(opt_omega)))]; 
H_desired(cut_freq)=[]; 

H_objective=(abs(H_desired)).‘2; 
[f]=simp1_3(H_objective,opt_omega(2)-opt_omega(1)) 


Of 9 Ae He He 6 2 ee He ee ee ee He AC Ee Ee ae AE HE HE He A HE eA EH A He Ke He Ae AC He He He ie 2 2 2 ie 2 ee 2 A 2 2 ae A 2 


J 
YJ Plot Generation 
lo ens arts 


% Saving a matrix of the optimization variables, objective function, and FRF in order to 


% plot how they changed during the optimization. 
J 
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if iter ==! 
del_vars = [del_kb del_kc del_cb del_cc]; 
H_desired_O = H_desired; 
del_f = [f]; 

else 
del_vars = [del_vars;del_kb del_kc del_cb del_cc]; 
del_H_desired = [H_desired_0; H_desired]; 


del_f = [del_f;f]; 
end 


clear hee H_desired_O H_desired 


fy 6 7 2 2 Ae 2 Ae 2 Ae 2 Ae 2 he 2 2 ee 2 he 2 he 2 he 2 Ae 2 ie 2 ie 2c 2c 2 2c 2 2 he 2c 2 2 2 2 ie 2 ie 2 ie 2 2c 2 ie 2 ie Ae 2 oe oe oe os 2 2 2 A KK 


%o 
% Calculation of plate and computer static defection 
D 


[stat_del_Kc,stat_del_Mc,stat_kb] = statdelta(cset,bset,kdof1,kdof2,1k,mdof,lm,kb0); 


[z_stat,Fred,Kred,Mred,Hstar0_dp] = statsynth(stat_hee,stat_kb,stat_omega ,... 
eset,iset,bset,stat_del_Kc,stat_del_Mc); 


Oy 75 2 As Ae 2 ee 2 2 Ae 2 Ae 2 Ae 2 he Ae ae he Ae 2 2 he Ac 2 fe 2c 2c 2 2 fe fe Ac 2 2 2c 2 he he 2 ae he 2 fe 2c 2 2c 2k fe 2c 2c 2k fe fe hc 2c 2 2c 2 2c 2c 2 2 oe 2s 2c 2 


Jo 
% Calculation of dynamic response due to shock 


[X_star,icnt] = ttmesynth(globalA,zset, Y, Ydot,del_Kc,del_Cc,del_Mc,kb,cb.... 
nstep,t_input(2),bset,cset); 


fprintf(‘icnt = %g',icnt) 
%odisp(‘Finish time synthesis’) 


Op 28 A Ae 2 Hee Hee He Ae He he 2 he 2 2 he 2 Ae he 2 ie 2 he 2 Ae 2 oie 2 ee ie ae 2 fe 2 fee ae he 2 he 2 he 2 ie 2 2 oe a eo KK A KK AK KK 


% 
% Calculation of normalized constraint values 


maxzp_sStat=-.08; maxzc_stat=-.03; 

maxkp_stretch=1; maxkc_stretch=1; 

maxzc_accel=30*386.4; 
Zp_stat=z_stat({1 3:6]); % pulls out the static displacements for the plate 
Zp_stat=-(max(abs(zp_stat))); % and determines where the maximum displacement 


%o 1S min command used because static displacement 
% 1s negative 


ZC_Stat=-(abs(z_stat(2)) - abs(z_stat(1))); % calculation of the relative static 


%o displacement between the plate and 
% computer 
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% calculation of the maximum relative dynamic displacements between the plate and the 
% base and the computer and the plate 


oO 
kp_stretch = X_star(2*nstep+1:6*nstep) - [Y;Y;Y;Y]; 
kp_stretch = max(abs(kp_stretch)); 
kc_stretch = X_star(nstep+1:2*nstep) - X_star(1:nstep); 
kc_stretch = max(abs(kc_stretch)); 


% calculation of the maximum absolute acceleration of the computer 


oO 
Xaccel_star = fddot(X_star((1:nstep)+(resp_index(1)-1)*nstep),t_input(2)); 
zc_accel = max(abs(Xaccel_star)); 


g(1)=zp_stat/maxzp_stat - 1; 
g(2)=zc_stat/maxzc_stat - 1; 
g(3)=kp_stretch/maxkp_stretch - 1; 
g(4)=kc_stretch/maxkc_stretch - 1; 
g(5)=ze_accel/maxzc_accel-1; 


save opt_data del_vars del_f iter del_H_desired g z_stat kp_stretch kc_stretch zc_accel ... 
X_Star Xaccel_star 
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